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SECTION  I 


INTRODUCTION 

The  area  of  Computational  Fluid  Dynamics  (CFD)  has  seen  tremendous 
progress  in  recent  years,  due  to  improvements  in  both  computational 
algorithms  and  computer  hardware.  The  improvements  now  make  it  possible 
to  obtain  steady  inviscid  solutions  for  very  complex  configurations  such 
as  modem  fighter  aircraft  and  the  external  weapons  used  on  these 
aircraft  (Reference  1).  However,  many  of  these  aircraft  and  weapons 
operate  in  flow  conditions  which  make  inviscid  solutions  an  inappropriate 
model.  Modern  fighter  aircraft  must  be  able  to  release  stores  and 
missiles  while  maneuvering  at  transonic  speeds.  The  analysis  of  these 
problems  requires  unsteady  calculations,  as  does  the  analysis  of 
aeroelastic  problems.  Many  modern  fighter  aircraft  have  internal  bays  in 
which  to  carry  their  weapons.  The  analysis  of  these  problems  requires 
that  a  highly  separated  unsteady  flow  field  be  solved.  Solutions  which 
include  the  effects  of  viscosity  are  necessary  to  appropriately  model 
these  flow  fields;  however,  for  these  types  of  configurations  viscous 
solutions  are  not  practical  with  current  CFD  codes,  due  to  the  large 
computer  run  times  and  large  memory  requirements.  The  emphasis  of  the 
work  presented  here  is  development  of  a  scheme  which  is  computationally 
efficient  and  accurate  enough  to  make  viscous  calculations  possible  for 
the  design  and  analysis  of  these  configurations  in  steady  and  unsteady 
flow  fields. 

The  approach  taken  is  an  extension  of  two  Drier  algorithms  de¬ 
scribed  by  Belk  :  and  Whitfield  ? .  The  schemes  used  in  this  work  are 
first,  a  flux  vector  split  scheme  based  on  Steger -Warming ^  splitting 
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and  secondly,  a  flux  difference  split  scheme  based  on  Roe’s  5 
approximate  Riemann  solver.  The  inclusion  of  the  viscous  terms  into 
these  Euler  algorithms  is  accomplished  as  described  by  Gatlin  s ,  where 
the  viscous  terms  have  been  added  to  the  algorithm  explicitly,  in  order 
to  save  the  cost  of  computing  and  storing  the  viscous  flux  jacobians.  The 
explicit  treatment  of  the  diffusive  terms  also  preserves  the 
computational  efficiency  peculiar  to  the  inviscid  algorithm. 

The  three-dimensional  Navier -Stokes  equations  in  differential  form 
are  presented  in  Section  II.  The  curvilinear  coordinate  transformation 
and  the  thin-layer  assumption  to  reduce  the  Navier-Stokes  equation  set 
are  also  presented. 

A  comparison  of  the  two  basic  implicit  algorithms  is  presented  in 
Section  III.  A  short  discussion  is  presented  which  compares  and 
contrasts  the  two  splitting  techniques,  flux  vector  splitting  (FVS)  and 
flux  difference  splitting  (FDS).  The  advantages  and  disadvantages  of  the 
two  algorithms  are  discussed,  with  the  discussion  based  on  the  one¬ 
dimensional  Euler  equations.  The  two  algorithms  are  then  extended  to  the 
three-dimensional  Euler  equations.  The  FVS  algorithm  has  been  presented 
in  some  detail  by  Whitfield  '  ( for  stationary  grids ) ,  and  Belk  -  ( for 
dynamic  grids)  and  therefore  is  only  summarized  here.  This  algorithm  has 
been  used  by  many  people  to  obtain  solutions  for  the  Euler  equations 
which  compare  well  with  experimental  data  for  both  steady  and  unsteady 

transonic  flow  problems.  Lijewski8  has  shown  excellent  steady  results 
for  a  complex  canard-body-tail  configuration  at  transonic  Mach  numbers 

and  moderate  range  of  incidence  angles,  while  Belk  9  has  shown 
insteady  results  for  oscillating  airfoils  and  wings.  However,  Gatlin  f 
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and  van  Leer111  have  shown  that  in  some  cases  the  FVS  scheme  is  too 
dissipative  to  give  good  numerical  results  for  viscous  calculations, 
especially  in  high  gradient  regions  such  as  boundary  layers.  Also 
presented  in  Section  III  is  a  discussion  of  the  FDS  scheme  based  on  Roe's 
method  of  evaluating  the  numerical  fluxes  for  the  three-dimensional  Euler 
equations.  Recent  works  by  van  Leer10  and  Chakravarthy  (References 
11,12,13,14)  and  many  others  have  discussed  implementations  of  this  FD6 
scheme.  However,  Whitfield,  et  al.3  and  Janus15  have  described  a 
hybrid  Roe  algorithm  which  has  been  used  to  obtain  excellent  inviscid 
results  for  complex  flow  fields16 . 

Section  IV  discusses  the  extension  of  the  inviscid  algorithms 
discussed  in  Section  III  to  the  thin- layer  Navier-Stokes  (TLNS) 
equations.  The  viscous  terms  are  added  to  the  explicit  side  of  the 
solution  equation  and  are  therefore  treated  as  source  terms  to  the 
inviscid  equations.  This  treatment  of  the  viscous  terms  reduces  the  cost 
of  normal  viscous  calculations  substantially  by  saving  the  cost 
associated  with  computing  and  storing  the  viscous  flux  jacobians. 
Additionally,  work  was  performed  to  allow  for  more  efficient  calculations 
of  unsteady  viscous  flows  using  the  time  lagged  viscous  terms.  One  of 
the  major  problems  with  unsteady  TLNS  calculations  is  the  severe  time 
step  size  limitations  which  must  be  imposed  to  maintain  a  reasonable 
Courant-Friedrichs-Lewy  (CFL)  number.  Efforts  to  execute  the  unsteady 
TLNS  code  for  oscillating  airfoils  and  wings  using  very  large  CFL  numbers 
are  discussed.  In  order  to  maintain  stability  at  these  large  CFL  numbers 
a  form  of  subiterations  similar  to  those  described  in  Reference  17  has 
been  Implemented.  The  result  of  subiterations  is  to  converge  the 
solution  at  a  given  time  step  before  preceding  to  the  next  step,  thus 
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eliminating  linearization  and  factorization  errors  and  increasing  the 
stability. 

Section  V  describes  the  steady  results  obtained  for  transonic  Mach 
numbers  and  moderate  range  of  incidence  angles  for  a  laminar  flat  plate, 
RAE  2822  airfoil  and  the  ONERA  M6  wing.  All  the  results  presented  in 
this  chapter  show  a  comparison  between  the  FVS  and  FDS  algorithms  for 
viscous  solutions  and  are  compared  with  experimental  data  where 
available. 

Section  VI  discusses  the  unsteady  results  for  both  oscillating 
airfoils  and  wings  in  transonic  flow  and  moderate  range  of  angles  of 
attack.  The  oscillating  airfoil  case  (NACA  0012)  is  compared  for  various 
time  steps  sizes  and  number  of  subiterations  while  the  oscillating  wing 
(Langley  Rectangular  Planform  Supercritical  Wing)  results  are  compared 

with  the  inviscid  results  of  Be lk  2  and  experimental  data. 

Conclusions  are  given  in  Section  VII. 
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SECTION  II 


NAVIER-STOKES  EQUATIONS 

2.1  Vector  Form  of  the  Navier -Stokes  Equations 
The  compressible,  three-dimensional  Navier -Stokes  equations  in 
Cartesian  coordinates,  written  in  strong  conservation  form  (without 
body  forces  or  external  heat  addition)  are 


■A  A  _ A  A 

.  3f  .  9g  ,  3h 

at  3x  3y  dz 


o 


(2. 


where  q,f,g,  and  h  are  dimensional  quantities  given  by 
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,  A  A  .  A  A  ~  A  ~  A  ~  A 

<e+p  )w~uT2x-VTzv_WT:.,+q^ 


and 


>\  1A,A4  /\-  v 

e  =  +  2P  u  +v  +w  ) 


The  first  row  in  Equation  (2.1 )  corresponds  to  the  continuity 
equation,  the  next  three  rows  are  referred  to  as  the  momentum  equations 
and  the  last  row  is  the  energy  equation.  In  the  strictest  sense  only  the 
momentum  equations  are  the  Navier -Stokes  equations;  however,  it  is  corrrron 
practice  to  refer  to  the  entire  set  of  five  equations  as  the  Navier- 
Stokes  equations. 

The  density  of  the  fluid  is  p,  the  pressure  is  p,  and  the  velocity 
conponents  in  the  Cartesian  coordinate  directions  x,y,  and  z  are  u,v. 
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and  ft  respectively.  The  relationship  between  &  and  £>  is  the  result 
of  the  perfect  gas  assumption  and  y  is  the  ratio  of  specific  heats. 


y 


CP 

rr 


(2.2) 


The  components  of  the  viscous  stress  tensor  are  given  as  7^. 

If  the  fluid  is  assumed  to  behave  as  a  "Newtonian  fluid",  that  is,  the 
stress  is  linearly  dependent  on  the  rates  of  strain  of  the  fluid  and  the 
coefficient  of  bulk  viscosity  is  negligible  so  that  the  relationship 

between  ^((coefficient  of  viscosity)  and  \( second  coefficient  of 
viscosity)  is 


(2.3) 


then  the  viscous  stress  tensor  in  Cartesian  coordinates  is  given  as. 


A  2  A  ,_A  A  A  . 

‘3  l'  <2«,  -v,  -V 


A 

Tyy  = 

A  2  A 


1  (l  <2Vv  -U,  -W.) 


^A  A  A 
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T„  *1  H  <2W,  -|1„  -Vy) 

A  A  A  ,  A  A 
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(2.4) 


M  '  1  A  y  A  /». 

T  =  T  =  11  (U  +  V  ) 

xy  yx  “  y  y 


A  A  A  ,  A  A  . 

-T„  -)1  <W„ 

A  A  A  ,  A  A 

T  =  T  =  11  (V  +  W  ) 

yz  zy  r  z  v 


and  the  heat-flux  vector  in  the  Cartesian  coordinate  directions  is  given 
as 


(2.5) 


7 


a  -  -  cf'l(  % 
"  ~ir~  iz 


where  Pr  is  the  Prandtl  number. 

It  should  be  noted  that  the  subscripts  on  the  velocity  ccrrponents 

and  on  the  terperature  T  denote  partial  differentiation, 
while  those  on  the  stress  and  heat  flux  terms  denote  Cartesian  coordinate 
directions . 

The  dimensional  quantities  (denoted  by  a  )  are  nondimensionalized 

using  frees tream  conditions  (denoted  by  »),  L  (the  reference  length 
used  in  the  Reynolds  number)  and. 
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(2.6) 


where 


A  A  A 

aOT  =  \JyRTh> 


is  the  frees tream  speed  of  sound. 

This  nondimensionalization  results  in  a  set  equations  of  exactly 
the  same  form  as  Equation  (2.1)  but  without  the  a  . 


+  af  +  +  aii  =  o 

3t  3x  3y  <3z 


(2.7) 
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2.2  Curvilinear  Coordinate  Transformation 


Equation  (2.7)  can  be  transformed  into  a  curvilinear  coordinate 
system  which  allows  the  use  of  body-conforming  grids  in  order  to  simplify 
the  implementation  of  boundary  conditions.  The  curvilinear  coordinates 
used  in  this  transformation  are 

£  =  £(x,y,z,t) 

T]  =  T](x,y,z,t)  (2.8) 

£  =  £(x,y,z,t) 

T  =  T^t) 

Note  the  use  of  time,  t,  to  define  the  curvilinear  coordinates.  This 
will  allow  for  the  time  dependent  motion  of  a  body.  The  details  of  a 
time  dependent  transformation  to  curvilinear  coordinates  cure  given  for 
the  Navier-Stokes  equations  in  Appendix  A. 


2.3  Thin-Layer  Navier-Stokes  Equations 

A  numerical  solution  of  the  Navier-Stokes  equations  can  require 
large  amounts  of  computer  time  and  memory.  One  assunption  often  made  to 
reduce  the  cost  of  a  solution  is  to  neglect  all  viscous  terms  which 
contain  derivatives  in  the  direction  parallel  to  the  body.  The  original 
development  of  this  idea  is  presented  in  Reference  18.  This  assunption 
is  justified  since  these  terms  are  substantially  smaller  than  the  terms 
with  derivatives  normal  to  the  body.  Also,  it  would  be  impractical  to 
think  a  grid  could  be  refined  enough  in  the  direction  parallel  to  the 
body  to  resolve  the  diffusive  terms  in  this  direction.  The  solution  of 
the  TLNS  set  of  equations  has  an  advantage  over  the  boundary- layer 
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equations,  which  assume  also  that  the  normal  pressure  gradient  is 
negligible.  Therefore,  the  thin  layer  Navier-Stokes  equations  (TLNS)  are 
capable  of  handling  flow  separation  and  reverse  flow  regions  with  no 
special  considerations  necessary. 

After  the  curvilinear  coordinate  transformation,  the  Navier-Stokes 
equations  have  the  form 

aQ  +  ^  +  9G  +  0H  =  o 
3t  3£  3 11  3C 

where  the  vectors  Q,F,G,  and  H  are  defined  in  Appendix  A,  Equation  (A.9). 
If  the  thin  layer  assumption  is  made  for  the  solid  surface  to  exist  on  a 
constant  p|  grid  line,  as  is  normal  for  a  "C  type"  numerical  grid,  then 
viscous  terms  inside  the  F  and  H  flux  vectors  cure  neglected  along  with  all 
the  terms  in  flux  vector  G  which  contain  derivatives  with  respect  to 
either  £  or  £ .  This  results  in  the  thin  layer  Navier-Stokes 
equations  as  follows 


3Q  +  M  +  3G  +  3H  _ 
3t  3£  3  ii  3C 


where. 
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and 
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Applying  the  thin  lay^-r  assumption,  the  viscous  stress  tensors 
described  in  Equation  (A. 10),  become 
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and  the  heat  flux  terms  become 
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SECTION  III 


DESCRIPTION  OF  SPLITTING  TECHNIQUES 


3 . 1  Basic  Irrplicit  Algorithm 

A  comparison  of  the  FVS  and  FIS  numerical  schemes  referenced  in 
Section  I  can  be  easily  understood  by  considering  the  one-dimensional 
Euler  equations 


3Q  +  olF 
,;)t  ,)J; 


0 


(3.1  ) 


An  implicit,  finite  volume  discretization  of  Equation  (3.1  )  may  be 
written  as 


AQ 11 
At 


=  o 


(3.2) 


where. 


AQ'1  =  Qn+1  -  Q 


and 


(F, 


I  *  I  /2 


1-1/ 


11+  I 


ah 


Since  in  the  curvilinear  coordinate  transformation  the  term  AH 
can  be  set  equal  to  one,  it  may  be  dropped  from  the  equation.  Equation 
(3.2),  when  solved  for  AQ  n  ,  becomes 

AQ 11  =  -At6pf"h  (3.3) 

Since  the  flux  F  is  a  non-linear  function  of  the  dependent  variables 
(p»pu,pv,pw,e)  then,  using  a  linearization  first  described  by 
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Briley  and  McDonald19  and  Beam  and  Warming  :o  where, 


Fn+I  =  Fn  +  A'^Q"  +  0(At: )  (3.4) 

the  first  order  accurate  scheme,  0(^t ),  based  on  the  backward  Euler 
differencing  formula.  Equation  (3.3)  is 

(I  +  At  )AQn  =  -At<%F"  (3.5a) 

while  a  second  order  accurate  scheme,  0(  /\t:‘ ) ,  based  on  the  three  point 
backward  differencing  formula  of  Equation  (3.3),  is 

(I  +  qAT  )^Q  "  -  -  ^j^F"  +  ^Q"~>  (3.5b) 

Note  that  the  dot(  )  indicates  the  5  operator  acts  on  the  entire  term 
( An^Q' ) ,  etc.,  and  the  term  A1',  is  defined  as  the  jacobian 
of  the  flux  vector 


A  = 


dF 

,|Q 


3 . 2  Flux -Vector  Splitting 

The  eigensystem  of  the  flux  jacobian  matrix  A  is  necessary  to 
achieve  the  desired  splitting.  The  flux  jacobian  matrix  can  be 
diagonalized  as  detailed  in  Reference  21 .  The  resulting  diagonal  matrix 
/\  has  the  eigenvalues  of  A  as  its  diagonal  elements, 

A  =  T,\T  1  (3.6) 

where  T  is  a  matrix  whose  columns  are  the  right  eigenvectors  of  A,  and 
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T  " 1 ' s  rows  are  the  left  eigenvectors . 


The  one-dimensional  Euler  equations  have  three  real  eigenvalues 
associated  with  the  jacobian  matrix  A.  The  eigenvalues  are  tne 
characteristic  velocities  of  the  system.  The  flux  vector  F  can  be 
written  as  the  sum  of  two  subvectors.  The  subvectors  are  written  with  a 
superscript  (+)  which  indicates  the  portion  of  the  flux  vector  resulting 
in  nonnegative  eigenvalues  and  superscript  (-)  which  indicates  the 
portion  of  the  flux  vector  resulting  in  the  nonpositive  eigenvalues 

F  =  F  *  +  F 

where 

which  allows  the  split  flux  jacobian  matrix  A'  to  be  defined 

A~  =  T/\  T  ■'  (3.9) 

Since  the  flux  vector  F  is  a  homogeneous  function  of  degree  one  in  Q 

then 


(3.7) 


(3.8) 


F 1  =  KQ  (3.10) 

The  splitting  of  the  flux  vector  F  described  above  is  the  flux  vector 
splitting  of  Steger  and  Warming  ■* . 

An  upwind  first  order  spatially  accurate  scheme  using  the  splitting 
described  can  be  constructed  as 
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(3.11  ) 


=  5  <F,  +  F,.,>  -  <iA|,„Q,.,  -  1*1,0,) 

where 

|  A|  =  A*  -  A' 

van  Leer10  has  pointed  out  that  many  numerical  fluxes  may  be 
written  in  the  form 

F1  +  i/2  =  ^(F,  +  Flt,)  -  1  D(Q,  ,Q1  +  |  )  *  ( Q,  * ,  -  Q,  )  (3.12) 

and  the  freedom  of  choosing  a  numerical  flux  formula  essentially  lies 
in  the  choice  of  the  matrix  coefficient,  D(Q,  ,Q1  +  |  )  of  the 
dissipation  term. 

The  dissipation  matrix  corresponding  to  the  Steger -Warming  splitting 
has  been  shown  to  contain  too  much  numerical  dissipation.  The  result  of 
this  additional  dissipation  is  a  smearing  of  shocks  and  contact 
discontinuities . 

3.3  Flux  Difference  Splitting 

The  flux  vector  split  scheme  described  above  has  been  described  by 

Barth22  as  being  a  "point -by-point  splitting"  or  in  other  words  a 
flux  vector  at  a  given  point  is  split  without  concern  of  the  state  of  its 
neighboring  points.  In  contrast  to  this  idea,  the  flux  difference  split 
scheme  referred  to  in  Section  I  is  based  on  the  "interaction"  of 
neighboring  points  through  an  approximate  solution  to  the  Riemann 
initial-value  problem.  The  FDS  scheme  used  in  this  work  was  based  on 
Roe's  approximate  Riemann  solver  which  uses  an  intermediate  state  of  the 
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dependent  variables  Q  to  determine  the  matrix  coefficient  A(Q )  of 
the  dissipation  term  described  in  Equation  (3.12).  The  numerical  flux 
formula  for  the  Roe  type  scheme  is 

f,.„2  -  2<f,  +  F,.|>  -  1*1,., -  Q.)  (3-'3> 

Roe  5  required  the  matrix  A(Q,  ,Q1  +  1  )  to  have  the  following  list  of 
properties,  which  he  calls  Property  U  since  it  ensures  uniform  validity 
across  discontinuities: 

(i)  It  constitutes  a  linear  mapping  from  the  vector  space  Q  to 
the  vector  space  F. 

RF 

(ii)  As  Q,  — >  Q(  +  1  -»  Q,  Then  A(Q,  ,Q|  +  1  )  —  A(Q),  where  A  =  ^ 

(iii)  For  any  Q, ,  Q1  +  |,  A(Qt,Q|  +  ))  •  (Q,  +  1  -  Q,  )  =  F(  +  1  -  F, 

(iv)  The  eigenvectors  of  A  are  linearly  independent. 

Prof^erty  (ii)  ensures  that  the  approximate  solution  approaches  the  exact 
solution  as  the  mesh  size  is  decreased.  Reference  10  points  out  that 
Property  (iii)  results  in  the  fact  that  if  the  eigenvalue  of  A(Q,  ,Q|  +  |  ) 
vanishes  then  the  corresponding  eigenvalues  of  the  dissipation  matrix 
|  A)  vanish,  since  the  eigenvalues  of  |A|  are  the  absolute  value  of 
the  eigenvalues  of  A.  Since  the  numerical  dissipation  vanishes,  this 
results  in  nice  sharp  modeling  of  shocks  or  contact  discontinuities  as 
opposed  to  the  Steger-Warming  FVS  scheme  in  which  the  eigenvalues  of  the 
dissipation  matrix  of  |  A|  become  discontinuous  whenever  the  corresponding 
eigenvalues  of  A  vanish. 

Let  \J,  j=1,2,...m  be  the  eigenvalues  of  the  Roe  matrix  A  and 
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rJ  be  the  right  eigenvectors  and  lf  be  the  left  eigenvectors  such  that 
1J  and  rJ  are  orthonormal  then 

|A  -  \'l|  =  0 

and 

(A  -  \'l)rJ  =  0 
lj(A  -  \JI)  =  0 

Jeffery23  (see  Chapter  2)  has  shown  that  the  change  across 
each  characteristic  curve  associated  with  a  specific  is  proportional 
to  the  right  eigenvector  r!  associated  with  that  eigenvalue  of  the  Roe 
matrix  A 


(3.14) 

(3.15a) 

(3.15b) 


Qi+i  Q] 


m 

CV.rJ 


(3.16) 


Jeffery  has  also  shown  that 


AdQ  1  =  \JdQ  > 


(3.17) 


Therefore,  extension  to  the  one-dimensional  Euler  pouations  results  in 


■ui  i 


F,  =  A(Q.,Qlt.)  *  (Q.  +  .  -  Q.) 


m 


=  r 


'I 


l-i 


m 

^df 1 

t-i 


(3.18a) 

(3.18b) 

(3.18c) 


where  df1  is  the  change  in  a  flux  vector  across  each  characteristic 

curve  associated  with  the  eigenvalue  v'  . 

A  numerical  formula  for  the  flux  vector  at  a  cell  interface  i+1/2  for 
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the  one-dimensional  Euler  equation  might  be 

m 

Fl  +  l/2  =  F,  +  <3.19) 

J-l  ' 

where  the  (-)  superscript  on  indicates  that  only  characteristic 
curves  associated  with  negative  eigenvalues  ( i.e. ,  nonpositive  slope  in 
Figure  3.1 )  are  used. 


t 


Another  flux  formula  is 


m 


F  =  F 
1+1/2  l+l 


U 


-  2>X'r 

i-i 


(3.20) 


where  the  (  +  )  superscript  on  'X  indicates  that  only  positive 
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eigenvalues  are  used.  Finally,  a  third  flux  formula  is  obtained  by 
averaging  Equations  (3.19)  and  (3.20) 

m 

f,.„2  =-2(f,  +f,m>  -i  Sa.IXV1  <3-21  ) 

J-l 

The  flux  formula  given  in  Equation  (3.21 )  can  be  interpreted  to  be 
in  the  form  described  by  Equation  (3.12)  by  letting 

$  =  |A|  (Q1  +  1  ~  Q,  )  (3.22) 

and  from  Equation  (3.16) 

m 

$  =  £|A|of  rJ  (3.23) 

J-i 


The  Roe  matrix  A  is  diagonalizable  just  as  the  flux  jacobian 
matrix  A  for  the  FVS  scheme  is  in  Equation  (3.6),  so  that 

|  A|  =  T(A+  -  A'  )T-'  (3.24) 

and 

m 

$  =  Vt(A+  -  A"  >T  (3.25) 

J-i  J 

Equation  (3.25)  can  be  shown  to  reduce  to  the  following  (see  Reference  3) 

m 

$  =  I.IXV’  (3.26) 

J-i 


Therefore  Equation  (3.21 )  can  be  written  as 


fi  +  i/2  ”  ^Fi 


+  fm>  ~ 


^1  Al 


1  +  1/2 


(Qi+,  -Q.) 
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(3.27) 


J 


The  extension  of  both  of  these  schemes  (FVS  and  FIDS)  to  higher  order 
accuracy  and  full  three-dimensional  equations  is  considered  next.  The 
discussion  will  deal  with  the  set  of  equations  as  though  they  were  purely 
hyperbolic  partial  differential  equations  (i.e.,  3-D  Euler  equations). 
This  is  possible  by  treating  the  diffusive  flux  terms  explicitly  and 
therefore  placing  them  on  the  right  hand  side  of  Equation  (2.9)  as  source 
terms.  The  details  of  this  explicit  treatment  of  the  viscous  flux  terms 
are  discussed  in  Section  IV. 

3.4  Flux-Vector  Splitting  for  Three  Dimensions 

Since  the  viscous  flux  terms  are  to  be  treated  explicitly  the 
remaining  terms  form  a  set  of  hyperbolic  partial  differential  equations 
characterized  by  a  limited  domain  of  dependence.  The  solution  at  a  point 
is  not  dependent  on  every  other  point  in  the  field  and  information 
travels  only  in  certain  characteristic  directions.  If  a  numerical  scheme 
is  selected  to  take  advantage  of  this  fact,  then  information  for 
numerical  differencing  is  taken  only  from  the  direction  from  which 
information  is  propagating  (i.e.,  upwinding).  The  use  of  upwind 
differencing  normally  alleviates  the  necessity  for  adding  numerical 
dissipation  to  smooth  the  results  near  discontinuities  such  as  shocks. 
This  becomes  apparent  when  one  considers  the  discussion  in  paragraph  3.1 
concerning  the  nature  of  the  dissipation  matrix  at  or  near 
discontinuities.  The  details  of  the  upwind,  finite  volume,  flux-vector 
split  scheme  for  the  three  dimensional  Euler  equations  have  been  reported 
many  times  (References  3, 7, 9, and  21),  and  thus  only  the  highlights  will 
be  provided  here. 

The  flux  jacobian  matrices  A,B,  and,  C  for  the  three-dimensional 
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unsteady  Euler  equations  are 


A  E 


rIF 

3Q 


The  eigenvalues  of  the  flux  jacobian  matrices  are  the  characteristic 
velocities  in  the  three  coordinate  directions.  The  flux  vectors  F,G,and, 
H  can  be  written  as  the  sum  of  two  subvectors  as  shown  in  Equation  (3.7) 
for  the  one-dimensional  equation.  The  subvectors  are  written  with  a 
superscript  (+)  which  indicates  the  portion  of  the  flux  vector  resulting 
in  nonnegative  eigenvalues  and  superscript  (-)  which  indicates  the 
portion  of  the  flux  vector  resulting  in  the  nonpositive  eigenvalues 


F  =  F  *  +  F‘ 

G  =  G  *  +  G  "  (3.28) 

H  =  H  +  +  H  “ 

Each  of  the  split  flux  vectors  are  linearized  the  same  as  described  in 
Equation  (3.4) 

F  +n+!  =  F  +n  +  A+AQn  +  0(&t2) 

F  'n+l  =  F  “n  +  A~AQn  +  0(At:  ) 

i  ,  -■  _  +n  +  l  _  _n*l  rT  *n*l  ,  _n*l 

and  similarly  for  G  ,G  ,H  ,  and  H 

The  split  flux  jacobians  are  now  defined  as  the  true  jacobians  matrices 
obtained  by  differentiating  the  positive  and  negative  flux  vectors 
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(3.29a) 

(3.29b) 


instead  of  the  splitting  defined  in  Equation  (3.9) 


A+ 


3F» 

9Q 


A" 


3F~ 

9Q 


(3.30) 


and  similarly  for  B  + ,  B  "  ,C  +  ,and  C  " .  The  derivation  of  these  split 
flux  jaoobian  matrices  are  given  in  excruciating  detail  for  dynamic  grids 
in  Reference  2. 

When  the  split  form  for  the  fluxes  given  by  Equation  (3.29)  is  used. 
Equation  (3.5),  when  expanded  to  three  dimensions,  becomes 

[i  +  s?a-+  6ntf+  6nB-+  6rc?+  5ccO]  AO1  =  -AtR1  (3.3ia) 

and 

R1  =  (5  e"  +  6  r  +  5  g"  +  5  gt  +  6  +  6rR~>  (3.31b) 

The  flux  vector  splitting  results  in  a  seemingly  simple  expression 
but  one  which  proves  impractical  to  solve  for  most  problems  of  any 
practical  size  due  to  the  wide  bandwidth  of  the  matrix  system.  Many 

different  factorizations  have  been  used  for  this  system  by  Whitfield  ' . 
The  factorization  used  in  this  work  was  the  two  factor  scheme,  in  which 
all  terms  corresponding  to  nonnegative  eigenvalues  (+)  are  lumped 
together  and  all  terms  corresponding  to  nonpositive  eigenvalues  (-)  are 
lumped  together. 

[i  +  At(S?a+-+  6rB+-+  6rcT-  )] 

[i  +  At(6?a"*  +  6rB"-+  6rc*  )]  Atf  =  -AtR1  (3.32) 
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The  factorization  error  created  is  no  larger  than  the  linearization  error 
already  inherent  in  the  scheme,  0(At2). 

The  solution  to  Equation  (3.32)  is  carried  out  as  follows 


[i  +  At  <5^.  +  6rB+-  +  5cc+-  )]  x1  =  -ATRn 

(3.33a) 

[i  +  +  8^8 '•  +  8rc-.  )]  x2  =  x1 

(3.33b) 

AQn  «  x2 

(3.33c) 

The  solution  of  Equation  (3.33)  results  in  two  passes  through  the 
ccnputational  field.  The  first  pass.  Equation  (3.33a),  solves  a  sparse 
block  lower  triangular  system  and  the  second  pass.  Equation  (3.33b), 
solves  a  sparse  block  upper  triangular  system.  This  results  in  a  very 
efficient  numerical  scheme  since  only  forward  and  backward  substitution 
are  required  and  not  a  matrix  inversion. 

Now  only  a  description  of  the  central  difference  operators  in 
Equation  (3.33)  is  needed  to  ccnplete  the  description  of  the  flux  vector 
split  scheme.  Since  the  scheme  is  finite  volume,  dependent  variables  are 
known  at  cell  centers;  however,  the  central  difference  operators 
represent  flux  differences  at  cell  faces.  The  value  of  the  dependent 
variables  at  cell  faces  are  determined  by  an  extrapolation.  For  a 
positive  extrapolation  a  value  at  the  cell  face  labeled  (i+1/2,j,k)  is 
determined  from 


+  o(Af) 


or 


Q\.w2  =  Q,  +  0<A£> 


(3.34a) 


(3.34b) 
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and  for  a  negative  extrapolation 

-&.> 

or 

Q‘,.,/2  '  c,.,  +  0<A£> 

The  first  order  extrapolations  are  used  for  the  difference  operators 
on  the  left-hand  (inplicit)  side  of  Equation  (3.33)  where  only  first 
order  spatial  accuracy  is  necessary  since  for  steady  results  — .0. 

The  second  order  extrapolations  are  used  on  the  right-hand  (explicit) 
side  of  Equation  (3.33)  to  give  a  second  order  spatial  accurate  scheme. 

In  order  to  evaluate  the  split  flux  difference  terms  in  the  residual 
vector.  Equation  (3.31b),  using  F  +  and  F  “  as  examples,  a  set  of  positive 
eigenvalues,  \(Q*),  and  positive  split  fluxes,  F,(Q+),  are  found  using 
the  dependent  variables  described  in  Equation  (3.34),  where  F^Q)  is  the 
component  of  the  flux,  F,  associated  with  the  eigenvalue  \  . 

Similiarly,  using  Equation  (3.35)  a  set  of  negative  eigenvalues,  \ ^ (Q“ ) , 
and  negative  split  fluxes,  F,  ( £r ) ,  are  determined. 

The  total  flux  through  a  cell  face  (i+l/2,j,k)  is  determined 
using 


(3.35a) 

(3.35b) 


F+  =  2  S  (q+  >  +  i\,(cr>i]  F,«r) 

i-i 

and 

F'  =  2  i  [^,(Q_)  -  IX,(Q')|]  f,(qt) 

l-l 


(3.36a) 


(3.36b) 
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The  numerical  scheme  described  here  has  been  labeled  BMULE  (Brown's 
Mule)  by  Whi+'f  ield  7 . 

3.5  Flux  Difference  Splitting  for  Three  Dimensions 

The  FDS  scheme  used  here  has  been  described  in  detail  by  Whitfield3 
for  application  to  the  three-dimensional  Euler  equations.  The  FDS  scheme 
of  Roe  is  based  on  an  approximate  Riemann  solver  which  solves  the  linear 
Riemann  problem  (shown  in  one-dimension,  where  A  is  the  jacobian 
matrix). 


—  +  A(Q,,Q,  ,  )  ^  =  0  (3.37) 

at  Wi  Vi+i  a* 

The  flux  formulas  discussed  in  paragraph  3.3  can  easily  be  extended 
to  three-dimensions  by  splitting  the  wave  into  the  three  curvilinear 
coordinate  directions  since  the  eigenvalues  are  dependent  on  the 
contravariant  velocities  normal  to  a  cell  face.  The  wave  motions  are 
split  based  on  these  eigenvalues  resulting  in  a  set  of  equations  in  which 
the  wave  motion  is  normal  to  the  cell  interface  just  as  though  they  were 
a  set  of  one-dimensional  equations. 

The  jacobian  matrix  A(Q(  ,Q1+(  )  is  evaluated  using  a  set  of 
special  averaged  values.  The  "Roe  averaged"  variables  for  the  three- 
dimensional  Euler  equations  are 


p = 

M  + 

u  ‘  4^  +  4^ 

.  Jp(v.  *  Iviv. 

Jp; +  JivC 

JiV”,  + 

w  = - —  — = — 

+  Jivi 


(3.38) 
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and 


(3.39a) 
(3.39b) 

where  h  is  the  total  enthalpy.  Roe  and  Pike  24  have  shown  that  these 
averaged  variables  are  the  only  ones  having  the  required  properties  to 
satisfy  the  Property  U  requirements.  The  details  of  the  derivation 
of  the  "Roe  averaged"  variables  are  given  in  Reference  5. 

An  implicit  algorithm  using  the  flux  difference  described  in 
Equations  (3.19)  and  (3.20)  is  obtained  by  linearization  of  the  first 
order  flux.  The  metrics,  denoted  by  M,  are  included  in  the  expansion  of 
the  flux  vector  and  the  jacobian  matrix  in  order  to  note  the  spatial 
location  that  the  metrics  are  to  be  evaluated.  Also,  noted  as  a 
superscript  is  the  time  step  to  which  the  metrics  correspond,  this 
becomes  important  for  calculations  using  dynamic  grids.  Equations  (3.19) 
and  (3.20)  are  written  in  the  following  form  using  the  same  development 
as  described  in  Reference  3,  where 

F n+l  =  Fn+1  +  A'  (Qn+I  -  Q n+l  )  (3.40a) 

1+1/2  1  1+1/2  l+l  I 

and 

F"+‘  =  Fn+1  +  A"  (Q n+l  -Qn+I)  (3.40b) 

1-1/2  I  i-1/2  I  1-1 

A  formal  linearization  of  Roe’s  numerical  flux  (Equation  (3.40)  is 

extranely  conplicated  as  recognized  by  Barth22 .  An  assumption  often 
made  to  simplify  the  linearization  is  to  neglect  the  spatial  derivatives 

of  the  Roe  matrix  A  .  If  this  assumption  is  made  then  the 


a2  =  -y  =  (y-i )  h  ~  ^(tr  +  v2  +  w2  > 
e  =  yzq  +  (u2  +  v2  +  w2 ) 
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linearized  flux  difference  at  cell  i  is 


Fn+1  -  Fn+1  =  Fn  -  Fn  +  A  AQn_A+  AQ"  + 

1+1/2  1-1/2  1+1/2  1-1/2  1+1/2  l+l  1-1/2  1-1 

[A(Q\  Mn+1  )  -  A"  -  A(Q n ,  M n+l  )  +  A+  ]£Q n  (3.41) 

1  1+1/2  1+1/2  1  1-1/2  1-1/2  U  I 


If  the  locations  of  the  metrics  in  the  jacobian  metrics  A  are 
neglected  then  Equation  (3.41)  results  in 


F  n+l 
1  +  1/2 


_-n+. 

1-1/2 


F  n 


-  F"  +  A'  A  Qn 

1+1/2  1-1/2  1+1/2  1+1 


A  +  A  Qn  + 

1-1/2  1-1 


[-  A 


1  +  1/2 


+  A  ]AQn 

1-1/2  1 


(3.42) 


or  if 

5(AAQ),  =  Altl/2  AQ^  -  A)  +  1/2  AQ" 

then 


-  n+l  -  n+l 

IT  _  C1 

1+1/2  1-1/2 


-  n  =-  n 

F  —  F 

i  +  1/2  1-1/2 


+  5a  ~AQn  +  5a  +AQn 


(3.43) 


(3.44) 


Hie  remaining  two  flux  terms,  G  and  H,  follow  in  a  similar  manner. 
Therefore,  Equation  (3.5),  when  expanded  to  the  three-dimensional  Euler 
equations  for  the  FD6  scheme,  results  in 


[I  +  At(5?A-+  5j;A‘+  5riB’+  5,^'+  5^C*+  5^C*)]  Acf  =  -At^  (3.45a) 
and 

F r  =  (5eF  +  5«F  +  5r)G  +  5r|G  +  5CH  +  5CH  )  (3.45b) 

Note  that  Equation  (3.45)  has  the  same  form  as  Equation  (3.31 )  from  the 
flux  vector  split  form  of  Equation  (3.5),  however,  the  flux  jacobian 
matrices  in  Equation  (3.45)  are  the  Roe  matrices,  and  the  jacobian 
matrices  in  Equation  (3.31 )  are  the  true  partials  of  the  positive  and 
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negative  flux  vectors.  Solutions  to  Equations  (3.31  )  have  been  obtained 
but  as  noted  in  Reference  3,  in  all  cases  tested  improved  results  and 
convergence  rates  were  obtained  from  solving  Equation  (3.31a)  but  using 
the  residuals  from  Equation  (3.45b) 

(3.46a) 

(3.46b) 

Therefore,  the  implicit  solution  from  the  FVS  scheme  is  used  with  the 
residuals  for  the  FD5  scheme. 

The  reason  for  the  improved  results  using  the  implicit  portion  of 
the  FVS  scheme.  Equation  (3.31a),  instead  of  the  implicit  portion  of  the 
FES  scheme.  Equation  (3.45a),  is  due  to  the  fact  that  approximated  flux 
jacobians  are  used  for  the  FES  scheme  while,  the  FVS  scheme  uses  the  true 
partials  of  the  flux  split  jacobians.  Any  approximations  for  these  flux 
jacobian  terms  have  nearly  always  degraded  the  convergence  rates  of  the 
solution. 

The  numerical  fluxes  discussed  in  paragraph  3.3  are  only  first  order 
accurate  in  space.  Higher  order  methods  used  here  (second  and  third 
order )  are  due  to  Osher  and  Chakra varthy  .  The  implementation  of  the 

higher  order  methods  is  described  by  Whitfield  3  (refer  to  Equation 
(54)). 

Using  Equation  (3.19)  for  the  first  order  flux,  a  second  or  third 
order  scheme  for  the  flux  at  cell  interface  i+1/2  may  be  obtained  by 


[I  +  At  <5^+  6n B*+  5nB-*+  orc*-+  srcO]  AQ1  =  -At*1 

and 

K1  =  (6?r  +  6  F  +  5rG+  +  5nG  +  6rH+  +  6  H  ) 
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adding  a  correction  term.  If  all  "Roe  variables"  and  metric  terms  used 
in  computing  the  eigenvalues  and  eigenvectors  are  evaluated  at  the  cell 
face  i+1/2,  then  the  higher  order  flux  formula  becomes 


**1  +  1/2  **^l  V1/2  + 


m 

V 

l-i 


0 


r  1 

|,l  +  1/2  1  +  1/2 


V  [L*(-1 ,1  )  -L" (3,1)] 


i-i 


H 


[L*(1,-1)  -  L~  ( 1 ,3)] }  r  1 

4  1  I  1+1/2 


(3.47) 


where 


and 


i  x  -J 

0  =  X  of,  ,  r 

J.l+p/2  1+1/2  f'~ 


<y 


j .  1-1/2 


J,  1  +  1/2 


11 ...  (Q,- Q,-, > 


1  +  1/2 


Of 


J. 1+3/2 


11  (Q,*rV 

1+1/2  11  1 

11  (Q1++-  Qi+. 

1+1/2  1  +  11 


(3.48) 

(3.49a) 

(3.49b) 

(3.49c) 


Two  flux  limiters  were  tested  in  this  work;  they  are  defined  as  the 
minted  and  superbee  limiters: 


L~  (l,n)  =  rninmod  ( 


O'  *  b 

j . 1+1/2 


O'  ) 

J. l+n/2 


minnod(x,y)  =  sign(x)  max  {0,min[|x|,y  sign(x)]} 

b  =  3J 

t  -  if. 


(3.50a) 

(3.50b) 

(3.50c) 


The  rninmod  limiter,  where  b  is  the  conpression  parameter,  uses  two 
arguments,  x  and  y.  When  these  arguments  are  of  opposite  sign  the  value 
returned  is  zero.  When  they  have  the  same  sign  the  value  returned  is  the 
smaller  absolute  value.  The  rninmod  limiter  uses  as  its  two  arguments  the 
unlimited  flux  value  at  a  cell  interface  and  compares  this  value  with  the 
product  of  the  compression  parameter  and  the  flux  value  at  the 


neighboring  interface  on  either  side.  When  these  values  have  the  same 
sign  than  the  value  returned  is  nearly  always  the  unlimited  flux  value 
since  the  compression  parameter  used  in  these  calculations  was  large  (two 
for  the  second  order  scheme  and  four  for  the  third  order  scheme).  Points 
in  the  solution  where  these  values  might  have  different  signs  are  the 
maxima  or  minima  of  fluxes  ( i.e. ,  shocks  or  strong  gradients).  In  these 
cases  the  unlimited  fluxes  can  easily  have  different  signs,  in  which  case 
the  minmod  limiter  returns  a  zero  and  the  solution  reverts  to  first 
order. 

The  superbee  limiter  used  in  conjunction  with  Equation  (3.47)  is 
simply 


(3.51a) 
(3.51b) 

(3.51c) 

where  (j  is  the  compression  parameter,  which  was  set  to  2  for  all 
results  shown  here. 

The  second  order  truncation  error  is  presented  for  several 
combinations  of  and  b  in  Table  I  from  Reference  13. 


L' (l,n)  =  I/  (n,l) 

J  ) 

L~(l,n)  =  anplim  (  gt  ,  Qt  ) 

1  J, 1+1/2  ] , i+n/2 

arplim(x,y)  =  sign(x)  {max  0,min[|x|,jjy  sign(x)]  , 
min[p|xj  ,y  sign(x)]} 


31 


SECTION  IV 


THIN-LAYER  NAVTER-STOKES  ALGORITHM 

4.1  Explicit  Treatment  of  Viscous  Terms 

The  diffusive  terms  may  be  written  as  separate  vectors  in  the  thin- 
layer  Navier-Stokes  equations.  Equation  (2.9)  becomes 


3Q  +  3 1  +  3(G-V  +  3«  = 

3t  3  c;  3ii  3C 


where 


(4.1  ) 


and  TTk  (k=x,y,or  z)  are  defined  in  Equation  (2.11  ). 

Since  the  Sv  vector  is  made  up  of  diffusive  terms  which  do  not 
have  che  characteristics  of  a  set  of  hyperbolic  partial  differential 
equations  in  which  a  set  of  characteristic  velocities  exist,  it  is 
impossible  to  split  the  diffusive  terms  into  subvectors  associated  with 
any  characteristic  velocities  and  impossible  to  upwind  difference  these 
terms.  Therefore,  the  inclusion  of  these  terms  into  the  left-hand  side 
of  either  Equation  (3.31)  or  Equation  (3.45)  would  destroy  the  efficiency 
of  the  upper  and  lower  triangular  matrix  solution,  by  requiring  a  central 
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differencing  of  these  terms.  In  an  effort  to  include  the  effects  of 
viscosity  and  yet  maintain  the  efficient  structure  of  the  algorithm 
Gatlirf  has  for  stationary  grids  neglected  the  inclusion  of  the  viscous 
flux  jacobians  in  the  left-hand  side  thereby  treating  the  viscous  terms 
explicitly.  The  result  is  the  viscous  terms  at  each  iteration  (n+1 )  are 
evaluated  using  the  information  (dependent  variables  and  metrics)  from 
the  previous  iteration  (n). 

An  implicit,  finite  volume  discretization  of  Equation  (4.1 )  with  the 
viscous  terms  becomes 

Aff  .  5F"*'  .  6G"*'  ^  SH"*'  .  f>Sv"*' 

&T  AC  All  AC  '  All 

However,  instead  of  the  linearization  of  Equation  (3.4)  for  the  viscous 
terms,  they  are  linearized  as 

Svn+I  =  Svn  +  O(At)  (4.3) 

The  remaining  fluxes  are  split  and  then  linearized  according  to  Equation 
(3.29)  and  the  equation  is  factored  into  two  factors  in  the  same  manner 
as  described  for  Equations  (3.33)  as 


[I  +  At(M+<  +  M+,+  bX- )] 

S  >1  l 

[i  +  At<5?a~-  +  6(b-*  +  5rcr-  )]  AO1  =  -At**1  (4.4) 

except  now  the  residual  term  includes  the  viscous  fluxes  which  contain 
derivatives  with  respect  to  the  normal  to  solid  surfaces. 
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(4.5) 


K1  =  <5^  +  h®  +  +  6r|Snv) 

The  residual  vector  in  Equation  (4.5)  can  be  used  in  either  a  FVS 
(Equation  (3.31b))  scheme  or  a  FDS  (Equation  (3.45b))  scheme.  This 
implementation  maintains  the  solution  algorithm  efficiency  of  the 
original  Euler  algorithm  (average  increase  in  computational  time  is  less 
than  10%)  while  including  the  viscous  effects  neglected  in  an  Euler 
algorithm.  Gatlin  6  has  shown  this  algorithm  to  give  reasonable 
engineering  solutions  for  high  Reynolds  number  flows  when  used  with  the 
FVS  scheme  for  stationary  grids.  Results  will  be  shown  in  this  work  for 
both  FVS  and  FDS  schemes  and  for  both  stationary  and  dynamic  grids.  The 
method  used  to  evaluate  the  viscous  flux  terms  has  been  described  in 
detail  by  Gatlin  5  and  is  not  repeated  here. 

4 . 2  Boundary  Conditions 

All  of  the  boundary  conditions  used  have  been  applied  explicitly. 

The  boundary  conditions  are  implemented  using  one  layer  of  phantom  points 
outside  of  the  computational  field,  which  results  in  a  first  order  in 
space  extrapolation  at  the  boundaries  and  enhances  the  vectorization  of 
the  computer  code.  The  phantom  points  at  far fields  are  set  to  enforce  a 
certain  condition  (supersonic  or  subsonic  inflow  or  outflow)  at  the  cell 
face,  which  is  on  the  boundary,  while  at  solid  surfaces  they  are  set  to 
enforce  the  no-slip  condition.  The  change  in  dependent  variable,  AQ1, 
is  set  equal  to  zero  for  all  boundaries  except  at  block  boundaries. 

All  far field  and  downstream  boundaries  used  characteristic  variable 
boundary  conditions  as  derived  in  Reference  21  for  stationary  grids  and 
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Reference  2  for  dynamic  grids.  Characteristic  variable  boundary 
conditions  are  designed  to  allow  information  to  flow  into  or  out  of  the 
computational  field  as  dictated  by  the  signs  of  the  eigenvalues. 


The  boundary  conditions  for  impermeable  surfaces  are  set  using  a  no¬ 
slip  implementation  of  the  zero  pressure  gradient  boundary  conditions. 
These  boundary  conditions  have  been  adapted  for  dynamic  grids  by  observing 

the  grid  speed  X  b  and  the  equations 


fp'P.d T?in  K  1  +  v 


(4.6a) 

(4.6b) 

(4.6c) 

+  w2  )  (4.6d) 


where,  the  subscripts  p,  f,  and  b,  denote  phantom  points,  field  points, 
and  boundary  points  respectively.  Figure  4.1  shows  graphically  the 
implementation  of  no-slip  condition  for  a  cell  centered  formulation. 


Figure  4.1  No-slip  boundary  conditions  for  dynamic  grids. 
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Finally,  the  block-block  boundary  conditions  were  derived  by  Be  Ik  2 
for  both  stationary  and  dynamic  grids,  lb  surrmarize  the  discussions  of 
different  block-block  boundary  conditions,  many  combinations  of  treatment 

for  AQn  and  dependent  variables  at  block  boundaries  exist.  Those 


used  for  the  results  obtained  hero  were,  two  phantom  points,  faQ 


approximated  (which  results  in  local  truncation  error  of 


t 


see  Belk  2 )  and  synchronized  dependent  variables. 


4 . 3  Turbulence  Modeling 

The  numerical  solution  of  the  Navier-Stokes  equations  ( Equation  2.1) 
for  turbulent  flows  require  the  viscosity  to  be  determined  by  the 
relation 


n  =  u0+nt  l4*7> 

where  n  is  the  molecular  viscosity  and  n  is  the  eddy  viscosity 
h  0  *1 

and  is  supplied  by  the  turbulence  model.  Similarly,  the  thermal 
conductivity  is  determined  by 

K  =  +  K,  (4.8) 

where  is  the  laminar  conductivity  and  is  obtained  from 

(4-9) 

while  n  is  the  turbulent  conductivity  determined  from 
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(4.10) 


where  Pt=0.72  and  Prl=0.90  for  air. 

Since  the  configurations  studied  here  are  airfoils  and  wings  with 
little  or  no  flow  separation,  the  Baldw in-Lomax  26  turbulence  model  was 
selected.  The  Baldwin-Lomax  turbulence  model  is  a  two- layer  algebraic 
eddy  viscosity  turbulence  model  which  uses  a  simple  algebraic  expression 
to  determine  a  value  for  m  in  the  inner  layer  nearer  the  wall 

*i,i 

and  n  in  the  outer  layer.  The  Baldwin-Lomax  model  has  the  advantage 

**  t  ,0 

of  not  requiring  the  determination  of  the  edge  of  the  boundary  layer  but 
relies  instead  on  the  vorticity  profile  to  determine  the  for  the 

inner  region: 


H,  |  =  Pl2lwl 

where 

1  =  KY[1  -  exp(-Y726)] 
and 

T  =  YJF^/P 

(jj  is  the  vorticity,  Y  is  the  normal  distance  from  the  solid  surface, 

«=0.4  is  the  von  Harman  constant,  r  is  the  maximum  of  r 

m  w 

and  t  ,  where  j  is  the  wall  shear  stress  and  t  is  the  maximum 

s  w  s 

shear  flow  in  a  local  velocity  profile.  The  results  of  the  selection  of  j 
in  the  flow  solution  is  discussed  by  Gatlin  6 . 


(4.11 ) 

(4.12) 

(4.13) 
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Fbr  the  outer  region  the  eddy  viscosity  is  given  by 


t  ,0 


«W.eb 


(4.14) 


where  Ccp=1.6  is  a  constant  and  K=0.0168  is  the  Clauser  constant. 
Fmax  is  the  maximum  of  the  profile  F(Y)  given  by 


and 


F  =  min  (Y  F  , 

w  'max  max1 


(4.15) 


F(Y)  =  Y|  oj |  [1  -  exp(-Yf/26)]  (4.16) 

is  the  y  value  at  which  Fm3X  occurs.  Q„k=0.25 
is  a  constant.  Fk|(,b  is  the  Klebanoff  intermittency  factor 
determined  from 


Fkleb(Y)  =  t1  +  5.5  (%tJ)6]-l 

hn  ax 


(4.17) 


where  C^.Ieb=0.3. 

Udlff  was  defined  by  Baldwin  and  Lcmax  as  the  difference 
between  the  maximum  and  minimum  velocity  magnitudes  in  a  profile 

U<j iff  =  Huy+  V2  +  W2)max  -  (4u:  +v 2  +  )mln 

4.4  Time  Accuracy  and  Subiterations 

Calculations  of  unsteady  aerodynamics  requires  a  time  accurate 
numerical  scheme  in  which  each  grid  cell  is  advanced  in  time  an  equal 
amount  (i.e.,  minimum  time  stepping).  The  size  of  a  time  step  is  limited 
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by  the  Gour ant -Friedrichs -Lewy  nunber  (CFL),  which  is  the  number  of  cells 
across  which  a  characteristic  wave  will  propagate  in  one  time  step.  If 
steady  state  solutions  are  of  interest  then  it  is  possible  to  update  each 
grid  cell  by  an  amount  in  time  which  places  that  cell  at  an  optimum  CFL 
number  for  convergence  (i.e. ,  local  time  stepping).  This  is  not  time 
accurate  since  each  grid  cell  will  be  at  different  time  levels,  but  this 
does  not  matter  when  a  steady  state  solution  is  being  sought  where  ACT— *0. 
All  the  steady  calculations  presented  here  resulted  from  local  time 
steps,  while  all  the  unsteady  calculations  resulted  from  minimum  time 
steps  in  which  each  grid  cell  was  updated  the  same  increment  in  time. 

The  use  of  minimum  time  steps  can  be  very  expensive  for  viscous 
solutions  since  grid  Is  near  solid  surfaces  are  very  small  compared  to 
grid  cells  in  outer  regions  of  the  flow  field  (6-7  orders  of  magnitude 
difference  in  cell  volumes).  This  means  very  small  time  steps  must  be 
used  over  the  entire  flow  field  for  minimum  time  stepping  to  maintain  a 
reasonable  maximum  CFL  in  the  grid  cells  near  solid  surfaces.  If  a 
maximum  CFL  of  100  is  occurring  in  the  grid  cells  near  a  solid  surface 
then,  most  of  the  flow  field  is  being  advanced  in  time  with  CFL 's  of  much 
less  than  one  (CFL  =  0.001-0.00001). 


Be lk  2  has  shown  good  inviscid  unsteady  results  using  this  FVS 
scheme  for  an  oscillating  airfoil  in  transonic  flow  using  a  maximum  CFL 
of  approximately  100.  To  obtain  a  maximum  CFL  of  100  on  the  inviscid 
grid  it  was  necessary  to  take  500  time  steps  per  cycle  of  oscillatory 
motion  of  the  airfoil.  To  obtain  the  same  approximate  maximum  CFL  (100) 
on  a  viscous  grid  would  require  approximately  500,000  time  steps  for  one 
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cycle  of  motion.  Obviously  unsteady  viscous  calculations  for  maximum  CFL 
numbers  on  the  order  of  100  are  very  expensive  and  are  not  practical  even 
on  today's  fastest  computers. 

If  unsteady  viscous  calculations  are  to  become  a  practical  tool  for 
use  in  design  and  analysis  of  transonic  fighter  aircraft  and  weapons  then 
it  will  be  necessary  to  reduce  the  number  of  time  steps  needed  for  a 
calculation  (i.e.,  increase  the  max  CFL  allowed).  With  this  implicit 
algorithm  a  linear  stability  analysis  indicates  an  unconditionally  stable 
scheme.  Gatlin  6  discusses  stability  concerning  the  explicit 
treatment  of  the  viscous  terms.  His  linear  stability  analysis  indicates 
the  scheme  maintains  its  unconditional  stability  for  high  Reynolds  number 
flows  only.  For  unsteady  problems  it  is  not  enough  to  simply  remain 
stable,  one  must  also  maintain  time  accuracy  when  using  these  large  CFL 
numbers. 

A  test  case  used  to  study  the  effect  of  CFL  on  stability  and  time 
accuracy  was  the  NACA  0012  airfoil  oscillating  in  pitch  ab^ut  its  quarter 
chord  point  in  transonic  conditions.  The  conditions  used  in  these 
calculations  were  taken  from  the  experiment  by  Landon27,  with  Mach 

number  of  0.755, "mean  angle  of  attack,  ,y  ,  of  0.016°,  an  unsteady 

“0  * 

angle  of  attack  amplitude,  ,  of  2.51°,  and  a  reduced  frequency 

'  U 

or  Strouhal  number,  k,  of  0.1628  based  on  chord  length, 

A  P* 

k  =  M_  (4.19) 

V, 

where 

(2)  =  frequency 


£  =  chord  length 
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and 


=  freestream  velocity 


The  angle  of  attack  of  the  airfoil  was  varied  according  to 

&(t)  =  o !  +  «u  sin(<2)t) 


(4.20) 


Unsteady  viscous  calculations  for  the  NACA  0012  airfoil  at  these 
conditions  were  made  for  different  time  step  sizes  to  determine  the 
effects  of  large  CFL  numbers  on  unsteady  viscous  calculations .  All  of 
the  NACA  0012  calculations  presented  were  obtained  using  a  ”C"  type 
numerical  grid  shown  in  Figure  1  with  221x40  points.  The  grid  was 
generated  using  the  Numerical  Grid  Generation  Code  -  EAGLE  written  by 
Thompson  (Reference  28).  A  great  deal  of  effort  was  expended  to  generate 
a  grid  which  had  as  much  orthogonality  as  possible  at  the  airfoil  surface 
(see  Figure  l.b),  since  the  FDS  scheme  was  found  to  be  very  sensitive  to 
orthogonality.  The  spacing  set  for  the  first  point  off  the  body  was 
0.000001  which  gave  approximate ly  15-20  grid  cells  inside  the  boundary 
layer  and  resulted  in  a  minimum  Y*  of  much  less  than  1  (actually  0.1 ) 
over  the  entire  airfoil. 


The  unsteady  calculations  were  performed  by  first  obtaining  a  steady 
solution  for  the  airfoil  at  the  mean  angle  of  attack,  (y  =0.016°. 

n 

The  motion  was  impulsively  started  with  the  angle  of  attack  varied 
according  to  Equation  (4.20).  The  entire  grid  was  oscillated  as  a  rigid 
body.  The  pressure  distributions  for  several  different  time  step  sizes 
are  presented  in  Figure  2.  These  pressure  distributions  are  "snap-shots" 
of  the  pressure  in  time  as  the  airfoil  is  in  motion.  The  airfoil  is  at 
the  6CP  of  motion  point  which  corresponds  to  an  increasing  angle  of 
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attack  through  2.20° .  The  small  time  step  size  DT  =  0.00102  was  used 
as  a  reference  time  step  to  compare  with  other  results  in  order  to  check 
for  time  accuracy.  As  the  time  step  increases  the  solution  deteriorates 
to  the  point  of  being  unusable  even  though  the  solution  did  not  "blow-up" 
and  continued  for  a  complete  cycle.  Table  II  suntnarizes  the  results  in 
Figure  2,  two  additional  time  step  sizes  have  been  included  in  Table  II 
for  comparison  purposes. 

The  results  shown  in  Figure  2  certainly  indicate  that  if  unsteady 
viscous  calculations  for  the  oscillating  airfoil  are  to  be  practical  or 
affordable,  additional  considerations  must  be  made  to  allow  calculations 
on  the  order  of  1000  time  steps  per  cycle  of  motion  or  less.  An  approach 
taken  to  help  in  this  cause  was  a  use  of  Newton  iterations  to  converge 
the  solution  at  each  time  step  before  updating  to  the  next  time  step. 
These  Newton  .iterations,  referred  to  here  as  subiterations,  have  been 
tried  successfully  by  others  (see  Reference  17),  but  not  with  the  two 
algorithms  described  in  paragraphs  3 . 1  and  3 . 2  and  not  with  the  explicit 
treatment  of  the  viscous  terms. 

The  implementation  of  subiterations  into  the  scheme  described  by 
Equations  (3.31  or  3.46)  is  a  simple  additional  term  added  to  the  right- 
hand  side  and  a  redefinition  of  A  CP  on  the  left-hand  side.  Using 
the  FVS  scheme  as  an  example  (the  implementation  is  exactly  the  same  for 
the  FDS  scheme),  one  has 


[I  +  £t(M+'  +  SnB**+  Sr0*'  ^ 

[I  +  At%A-+  5  BT-+  !)rC-  )]  A01'  =  ~Atr1'  (4.21a) 
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-At  =  -AtC2^-  +  <5^  +  5ncr  +  8^’  +  6^)]  (4. 21b) 

where  ^Gp  is  defined  as 

AQf  =  Q^1  -  Qf'  (4.21c) 

where  p  is  the  subiteration  count.  When  p=0,  Q  p  =  Qn,  and  the 
system  reverts  to  the  noniterative  scheme  of  Equation  (3.31 ).  At 

convergence,  Q  p  — *  Q  n+1 .  Note  that  in  this  case  the  viscous  terms 
are  not  time  lagged.  The  normal  operating  procedure  was  to 
simply  set  the  number  of  subiterations  desired  rather  than  check  the 
convergence  level  after  each  subiteration. 

Qretga  and  Rheinboldt29  write  the  Newton  iteration  in  the  form 

XT1  =  yP  -  )''F(X)P  p=0 ,1,2,...  (4.22) 

where  k(p)  is  an  integer  less  than  or  equal  to  p.  When  k(p)  is  less  than 
p  then  the  F'(X)  is  re-evaluated  less  than  each  iteration  which 
translates  into  "freezing"  the  flux  jacobians.  When  k(p)  =  p  the 
jacobians  are  evaluated  each  iteration  and  Equation  (4.22)  results  in  the 
normal  Newton  iteration 

XT'  =  X?  -  F,(Xp)_lF(X)p  p=0, 1,2,...  (4.22) 

When  k(p)  =  0  the  jacobians  are  never  updated  and  Equation  (4.22)  is  then 
referred  to  as  the  simplified  Newton  method 
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XP+I  =  -  F'(X°  )  'F(X)P  p=0, 1,2,  —  (4.23) 

The  significance  of  freezing  the  jacobians  and  the  resulting 
sinplified  Newton  method  is  a  tremendous  savings  in  corrputations  since 
the  cost  to  recalculate  the  flux  jacobians  is  between  40%  -  45%  of  the 
cost  of  doing  an  entire  iteration.  Both  methods  were  tried  with  no 
noticeable  differences  in  the  results,  therefore,  only  the  sinplified 
Newton  method  results  will  be  reported  on. 

Pressure  distributions  for  the  NACA  0012  at  60P  motion  cure  shown 
in  Figure  3  for  two  time  step  sizes.  The  smaller  time  step  size,  DT  = 
0.0511,  ( see  Table  III  for  a  description  of  the  conditions  for  the 
subiteration  check  cases)  shows  some  moderate  wiggles  for  the  zero 
subiterations  while  for  only  four  subiterations  the  solution  is  much 
smoother  and  is  suitable  for  determining  the  lift  and  moment  coefficients 
for  use  in  design  and  analysis  studies.  However,  for  the  larger  time 
step  size,  OT  =  0.2044,  the  original  solution  without  subiterations  is 
very  irregular  and  even  after  32  subiterations  the  solution  still  shows 
(Figure  4)  a  significant  number  of  wiggles.  The  solution  could  be 
smoothed  but  would  require  numerous  subiterations  due  to  the  slow 
convergence  rates  for  the  subiterations  and  would  not  be  cost  effective. 

The  idea  of  using  subiterations  for  unsteady  calculations  was  to 
converge  the  solution  at  each  time  step  and  therefore  allow  for  much 

larger  time  steps  sizes.  However,  Swanson  and  Turkel30  have 
reported  on  the  problems  of  converging  solutions  with  numerical  grids 
designed  for  viscous  flows  which  have  a  large  variation  in  grid  cell 
aspect  ratio  between  grid  cells  near  a  solid  surface  and  grid  cells  in 
the  far  field.  A  solution  for  the  convergence  rate  problem  in 
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subiterations  which  worked  well  was  to  use  the  idea  for  accelerating 
steady  state  convergence,  that  is  local  time  stepping. 

Local  time  steps  used  for  the  subiterations  seemed  to  have  a 
tremendous  inpact  on  the  convergence  rate  for  the  subiterations.  The 
local  time  steps  were  implemented  such  that  the  overall  time  accuracy  of 
the  minimum  time  step  iterations  was  preserved.  Equation  (4.21b)  then 
becomes 

-At  -  -AT,oca, +  <6ff  ♦  &„<?  +  +  5,=t>]  (4.25) 

*  min 

The  result  of  using  local  subiterations,  Equation  (4.25),  was  a 
significant  increase  in  the  convergence  rate  for  subiterations  as  long  as 
the  CFL  for  the  local  subiterations  was  small  enough  to  maintain 
stability.  Numerical  experiments  showed  that  the  local  CFL  must  usually 
be  less  than  one.  This  fact  would  eventually  become  the  Achilles  heel 
for  local  subiterations  since,  determining  the  optimum  local  CFL  number 
which  allowed  the  Newton  subiterations  to  remain  stable  and  still 
converge  rapidly,  proved  to  be  extremely  difficult.  A  method  guaranteed 
of  selecting  the  optimum  local  CFL  number  was  not  found;  however, 
provided  a  proper  value  for  the  local  CFL  was  found  the  local 
subiterations  converged  faster  than  the  minimum  time  steps. 

The  pressure  distributions  in  Figure  5  show  very  smooth  results  for 
tine  step  size,  CT  =  0.05111,  using  local  time  stepping  subiterations 
even  for  only  two  subiterations.  The  results  in  Figure  6  show  much 
smoother  pressure  distributions  for  time  step  size,  DT  =  0.2044,  than 
did  the  minimum  time  steps  ( Figure  4 ) .  Even  though  there  is  a 
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significant  amount  of  ringing  near  the  shock,  this  solution  would  seem  to 
be  sufficiently  smooth  to  give  reasonable  engineering  answers  for  lift 
and  moment  coef ficients. 

Figure  7  shows  a  comparison  for  the  time  step  size,  DT  =  0.05111, 
with  4  minimum  subiterations  and  4  local  subiterations.  The  results  are 
very  similar  for  the  two  cases.  However,  both  subiteration  cases  show 
the  shock  wave  farther  aft  than  does  the  smaller  time  step  size  with  no 
subiterations.  Since  the  shock  is  still  building  in  strength  and  moving 
aft,  the  larger  time  step  size  shock  position  is  said  to  be  "leading"  the 
smaller  time  step  size.  This  situation  was  consistently  observed  for 
these  large  time  step  size  calculations. 

The  time  accurate  calculation  of  unsteady  aerodynamics  would  at 
first  glance  seem  to  require  a  second  order  temporal  accurate  scheme. 

However,  Be lk 2  has  shown  that  for  Euler  calculations  on  oscillating 
airfoils  and  wings  there  is  very  little  difference  between  solutions 
obtained  with  a  second  order  accurate  scheme  and  those  obtained  using  a 
first  order  accurate  scheme.  The  extension  of  the  Euler  algorithm  to 
include  the  diffusive  terms  (paragraph  4.1 )  requires  the  diffusive  terms 
to  be  time  lagged  and  results  in  a  first  order  scheme  regardless  of  the 
time  discretization  used.  For  example,  consider  the  linearization  as 
described  in  Equation  (3.4) 

r"+l  =  F1  +  AnA£P  +  0(^t2)  (3.4) 

When  this  linearization  is  used  in  conjunction  with  the  second  order 
time  discretization  given  in  Equation  (3.5b)  the  result  is  a  second  order 
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accurate  scheme.  When  Equation  (3.4)  is  used  in  conjunction  with  the 
first  order  time  discretization  given  in  Equation  (3.5a),  the  result  is  a 
first  order  accurate  scheme.  When  the  viscous  terms  are  added  explicitly 
using  the  linearization  given  in  Equation  (4.34),  the  formal  time 
accuracy  is  dropped  to  0(^t)  regardless  of  the  order  of  the  time 
discretization  used,  either  Equation  (3.5a)  or  (3.5b). 

The  use  of  subiterations  would,  it  appears,  alleviate  this  problem 
of  time  accuracy,  since  at  convergence,  the  linearization  and 
factorization  errors  go  to  zero.  Therefore,  using  the  second  order  time 
discretization  described  in  Equation  (3.5b)  with  subitex.a Lions,  gives  a 
second  order  accurate  in  time  scheme  at  convergence.  A  comparison  of  a 
first  and  second  order  differencing  with,  DT  =  0.0511,  is  shown  in  Figure 
8  using  16  minimum  subiteration  time  stepjs,  which  resulted  in  an  order 
of  magnitude  reduction  in  the  residuals.  There  are  seme  differences, 
with  the  major  difference  being  the  shock  location.  Since  the  shock  at 
this  point  in  the  motion  is  still  strengthening  and  moving  aft  very 
rapidly,  then,  as  would  be  expected,  the  first  order  shock  position  lags 
the  second  order  shock  position. 

The  result  of  these  investigations  was  to  help  select  the  time  step 
sizes,  number  and  type  of  subiterations  to  be  used  for  calculations  of 
the  oscillating  airfoil  and  wings  for  several  cycles  of  motion.  Minimum 
time  step  subiterations  were  used  because  of  the  problems  of  selecting 
the  local  CFL  for  the  local  time  step  subiterations.  A  time  step  size 
considered  appropriate  for  inviscid  calculations  was  used  with  4-5 
subiterations,  which  resulted  in  a  good  tradeoff  between  cost  and 
accuracy.  These  results  will  be  presented  in  Section  VI,  but  first  a 
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comparison  of  the  EVS  and  FDS  schemes  for  steady  viscous  calculations 
will  be  presented  in  Section  V. 
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SECTION  V 


STEADY  RESULTS 

The  description  of  the  splitting  techniques  in  Section  III  indicated 
the  FDS  scheme  based  on  Roe's  approximate  Riemann  solver  should  out 
perform  the  FVS  scheme  (BMULE)  based  on  the  Steger -Warming  splitting. 

Vein  Leer  and  Gatlin  have  pointed  out  that  the  numerical  dissipation  for 
the  TVS  scheme  is  too  large  to  allow  the  accurate  modeling  of  flow 
discontinuities  such  as  shocks  and  boundary  layers.  In  order,  to 
document  the  effects  of  numerical  dissipation  for  viscous  flow  solutions, 
several  steady  configurations  were  tested  to  compare  the  FDS  scheme  with 
the  FVS  scheme. 

5.1  Flat  Plate  Boundary  Layer 

The  first  configuration  was  a  laminar  flat  plate  boundary  layer 
calculation  for  Mach  =  0.5  and  a  Reynolds  number  based  on  plate  length  of 
Re  =  10,000.  Flow  solutions  for  a  laminar  boundary  layer  on  a  flat  plate 
were  obtained  using  two  grids.  Both  grids  used  70  points  in  the 
freestream  direction  with  1 0  points  in  front  of  the  plate  stretched  from 
the  outer  boundary,  located  4.0  plate  lengths  in  front  of  the  plate  to 
the  plate  leading  edge.  The  next  50  points  were  uniformly  distributed  on 
the  plate,  and  the  last  10  points  were  placed  aft  of  the  plate,  stretched 
from  the  plate  trailing  edge  to  5  plate  lengths  downstream  of  the  plate. 
The  two  grids  differed  in  the  number  of  grid  lines  in  the  normal 
direction  frem  the  plate  and  in  the  spacing  set  for  the  first  point  off 
the  plate.  The  first  grid  used  20  points  in  the  normal  direction  with 
the  first  point  off  the  plate  set  at  0.002  plate  lengths  (fine  grid), 
while  the  second  grid  used  15  points  in  the  normal  direction  and  the 
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first  spacing  was  0.02  plate  lengths  (course  grid). 


The  laminar  boundary  layer  profiles  are  presented  in  Figure  9  for 
the  Roe  scheme  on  the  course  grid,  and  EMJLE  on  both  the  course  and  fine 
grids.  The  calculations  are  shown  compared  with  the  Blasius  solution  for 
a  laminar  flat  plate.  The  Roe  scheme  matches  very  well  with  the  Blasius 
solution  using  the  course  grid.  Note  that  the  Roe  scheme  has 
successfully  modeled  the  laminar  boundary  layer  profile  with  only  three 
mesh  points  internal  to  the  boundary  layer.  The  fine  grid  results  with 
the  BMULE  scheme  do  not  compare  as  well  as  the  Roe  scheme.  The 
additional  grid  refinements  did  not  improve  the  BMULE  results.  The  BMULE 
scheme  has  too  much  numerical  dissipation  to  capture  the  laminar  boundary 
layer  profile,  even  on  the  fine  grid. 

Calculations  with  the  two  schemes  were  also  completed  for  a 
turbulent  flat  plate  boundary  layer.  While  the  Roe  scheme  still  out 
performed  the  BMULE  scheme,  the  differences  were  not  so  dramatic.  The 
Roe  scheme  was  able  to  model  the  turbulent  boundary  layer  with  a 
relatively  course  grid  once  again.  However,  it  was  necessary  to  ensure 
that  at  lease  one  grid  point  was  inside  the  laminar  sublayer  of  the 
turbulent  boundary  layer.  Therefore,  the  Roe  scheme  or  the  BMULE  scheme 
requires  a  spacing  for  the  first  grid  point  off  the  body  to  be  such  that 
a  minimum  Y*  ( Equation  4.13)  of  less  than  2-3  occurs .  This  criterion 
was  used  in  the  remaining  calculations  to  ensure  the  turbulent  boundary 
layers  were  being  satisfactorily  modeled. 

5.2  RAE  2822  Airfoil 

The  RAE  2822  airfoil  was  used  as  a  check  case  for  the  thin- layer 
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Navier -Stokes  algorithms  using  both  Roe  and  BMCJLE.  The  RAE  was  selected 
because  of  the  experimental  data31  which  exists  for  pressure,  skin 
friction,  and  velocity  profiles. 

The  RAE  2822  is  a  supercritical  airfoil  as  shown  in  Figure  10,  with 
a  design  Mach  number  of  0.66.  A  221  x  31  x  2  grid  was  used  to  compute 
the  flow  for  the  case  presented.  The  31  grid  points  were  stretched  from 

1  x  10'6  spacing  for  the  first  point  off  the  body  (close  enough  to 
guarantee  a  minimum  Y*  <  1 ),  5  chord  lengths  to  the  front  far field 
boundary,  and  10  chord  lengths  to  the  top  and  bottom  far  field  boundaries. 
The  downstream  boundary  was  set  at  15  chord  lengths  with  30  points  in  the 
wake.  The  remaining  points  (160)  were  spread  around  the  airfoil,  with 
clustering  at  the  leading  and  trailing  edges. 

Flow  solutions  were  obtained  for  the  RAE  using  the  grid  described 
above  for  four  algorithm  combinations.  The  Roe  scheme  was  used  with 
mirmod  limiter  in  both  the  second  and  third  order  spatially  accurate 
mode,  and  with  the  superbee  limiter,  which  is  second  order  accurate. 

The  BMULE  scheme  was  used  to  compare  with  each  of  the  Roe  scheme 
algorithms.  The  flow  conditions  for  the  computation  were  taken  at  Mach  = 
0.73  and  angle  of  attack  =  3.19°.  The  flow  conditions  were  corrected 
to  account  for  flow  angularity  and  wall  effects  to  the  condition 

recommended  by  Cook  et  al31.  The  calculations  were  at  Mach  =  0.734 
and  angle  of  attack  =  2.79 p.  The  Reynolds  number  for  both  the 

experiment  and  the  calculations  was  6.5  x  1 0  6 ,  based  on  chord 
length. 
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Figures  1 1  .a  -  1 1  -d  show  the  ccnparison  for  the  Roe  mimed  limiter 
with  2nd  order  spatial  accuracy  and  the  EMCJLE  2nd  order  accurate  scheme. 
Figure  1 1 .a  shows  the  pressure  distribution,  while  11  .b,  and  1 1 .c  show 
conparison  of  velocity  profiles  at  two  stations  on  the  upper  surface  of 
the  airfoil.  Figure  ll.b  shows  the  velocity  profile  at  32%  of  the  chord, 
which  is  well  ahead  of  the  upper  surface  shock,  while  Figure  1 1 .c  shows 
the  velocity  profile  at  57%  chord,  which  is  very  near  the  upper  surface 
shock.  Figure  1 1 .d  shows  a  conparison  of  skin  friction  between  the  two 
algorithms  and  experimental  data  for  the  upper  surface. 

Figures  12. a  -  12. d  and  Figures  13. a  -  13. d  follow  the  same  pattern 
as  Figures  11. a  -  1 1  .d,  with  Figures  12  comparing  the  Roe  minmod  3rd 
order  scheme  with  the  BMULE,  and  Figure  13  comparing  the  Roe  super bee  2nd 
order  scheme  with  the  BMJLE  scheme. 

Each  of  the  algorithm  cases  compare  very  favorably  with  the  pressure 
data  t Figures  (11  -  1 3 1 . a ) .  Roe  schemes  match  better  in  the  shock  regions 
due  to  the  numerical  dissipation  in  the  EMULE  scheme  smearing  the  shocks. 
The  Roe  superbee  scheme  was  the  only  one  to  accurately  model  the  leading 
edge  expansion  for  the  RAE.  However,  the  Roe  super bee  did  not  do  as  well 
in  the  shock  region  due  to  superbee  predicting  a  small  shock  induced 
boundary  layer  separation  bubble,  which  was  not  evident  in  the 
experimental  data  or  in  the  other  calculations. 

The  velocity  profiles  shown  in  Figures  ( 1 1 , 12,13 )b-c  are 
plotted  as  UC/Ue  (contra variant  velocity/contra variant  velocity  at  the 
edge)  versus  the  normal  coordinate  distance  from  the  upper  surface.  The 
edge  of  the  boundary  layer  was  determined  to  be  the  cell  at  which  the 
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maximum  If  and  contravariant  velocity  values  occurred.  When  these 
maximum  values  did  not  occur  in  the  same  grit'  cell,  a  simple  averaging 
was  used  to  determine  Ue.  As  is  evident  in  all  the  plots,  the  Roe  scheme 
is  far  superior  to  the  EMULE  scheme.  The  BMULE  scheme  consistently 
predicts  a  thicker  boundary  layer  than  that  shown  in  the  experimental 
data.  This  again  is  attributed  to  the  numerical  dissipation  in  the  BMULE 
scheme.  The  superbee  limiter  as  discussed  above  predicts  some  flow 
separation  behind  the  shock  which  is  made  more  evident  by  observing  the 
velocity  profile  for  57%  chord  position.  Figure  13.c.  In  all  other  cases 
the  Roe  scheme  compares  very  well  with  the  experimental  velocity 
profiles,  and  this  is  true  for  these  calculations  even  though  the  grid 
contained  10  less  points  in  the  normal  direction  than  did  Gatlin's 
6  ,  which  would  be  considered  by  some  as  already  being  a  sparse  grid. 

(See  References  32  and  33.) 

A  comparison  of  skin  friction  between  the  Roe  scheme,  EMULE  scheme, 
and  experiment  are  shown  in  Figures  (11-13)d.  Both  the  Roe  scheme  and 
BMULE  tend  to  overpredict  the  skin  friction  when  compared  to  experimental 
data.  Most  of  the  calculations  match  well  with  the  experimental  data  up 
to  the  shock  location,  but  both  schemes  tend  to  greatly  overpredict  the 
skin  friction  behind  the  shock.  The  skin  friction  data  near  the  trailing 
edge  for  Figures  ll.d  and  12.d  was  lost  due  to  a  data  file  transfer 
error.  Once  again  the  evidence  of  superbee's  separated  region  behind  the 
shock  is  evident  in  Figure  13.d. 

Even  though  the  Roe  scheme  gives  improved  results  in  almost  all 
cases  shown,  the  improvements  do  not  come  cheaply.  Figure  14  shows  a 
convergence  history  for  the  two  schemas.  Obviously,  the  BMULE  scheme 
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converges  much  faster  than  does  the  Roe  scheme,  so  much  so  that  after 
only  800  iterations  at  a  CFL  =  15,  the  BMULE  residuals  are  nearly  an 
order  of  magnitude  smaller  than  are  the  Roe  scheme  residuals.  Numerical 
experimentation  showed  that  800  iterations  was  the  minimum  necessary  to 
obtain  engineering  answers  (after  which  the  solution  did  not  change)  for 
the  Roe  scheme  while  only  500  was  necessary  for  BMULE.  This  translates 
into  many  fewer  iterations  necessary  to  arrive  at  the  same  convergence 
level  for  the  BMULE  scheme.  When  this  is  added  to  the  fact  that  the  Roe 
scheme  has  a  higher  operations  cost  than  does  BMULE,  nearly  20%,  since 

BMULE  runs  at  6.185  x  10'5  CPU  sec/iteration/point  while  Roe  runs  at 

7.417  x  10-5  CPU  sec/iteration/point,  it  becomes  evident  the  Roe  scheme 
is  much  more  expensive  to  use. 

5.3  ONERA  M6  Wing 

The  ONERA  M6  is  a  syrrmetric  12%  thick  airfoil  section  with  a  sweep 
angle  of  30  degrees.  The  wing  is  tapered  with  a  taper  ratio  of  0.56  and 
has  an  aspect  ratio  of  3.8.  Extensive  wind  tunnel  test  data  exist  for 
the  ONERA,  in  particular  pressure  data  for  transonic  flow  conditions 
(Reference  34).  The  ONERA  is  used  here  to  compare  the  BMULE  scheme  and 
the  Roe  scheme  with  experimental  data  for  a  steady  three-dimensional 
configuration  at  Mach  =  0.84,  angle  of  attack  =3.06  degrees,  and  Re  = 

2.6  x  10  6. 

The  grid  used  in  these  calculations  was  a  111  x  40  x  25  "C-O"  type 
grid.  The  grid  was  generated  using  a  distribution  similar  to  that  used 
for  the  NACA  0012  airfoil  grid  described  in  paragraph  4.4.  The  upper  and 
lower  airfoil  sections  in  this  grid  were  generated  independently  with  the 
first  point  in  the  grid  (1,1,1)  being  at  the  wing  leading  edge.  The  grid 
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used  a  spacing  1x10s  for  the  first  spacing  off  the  body,  which 
resulted  in  a  min  Y*  <  1  over  the  entire  grid.  The  outer  boundaries 
were  extended  to  10  chord  lengths  in  all  directions.  Figure  15  shows  a 
closeup  of  two  K  planes  (constant  spanwise  location)  and  the  airfoil 
section.  The  grid  was  generated  by  letting  the  K=1  plane  be  the  plane  of 
grid  points  on  the  upper  wing  surface  at  the  root  section.  The  K-planes 
were  then  distributed  linearly  along  the  upper  surface  to  K=8  at  the 
upper  surface  tip  (Figure  I5.b).  Planes  9  through  17  were  then  rotated 
in  a  circular  arc  to  model  the  wing  tip,  while  K-planes  18  through  25  on 
the  lower  surface  were  distributed  linearly  from  the  tip  to  the  root. 

This  wraparound  wing  tip  grid  shown  in  Figure  15.c,  allows  the  modeling 
of  the  wing  tip  as  it  existed  in  the  wind  tunnel  model. 

The  test  case  was  run  for  four  algorithms  just  as  for  the  RAE 
airfoil.  The  BMULE  algorithm  was  compared  to  the  Roe  algorithm  using  the 
2nd  order  mirmod  limiter,  3rd  order  minmod,  and  2nd  order  superbee 
limiter.  The  test  case  was  a  transonic  condition  which  results  in  a 
double  shock  configuration,  which  is  evident  in  Figure  16. a  at  near 
midspan  point  (44%).  Figure  16.b  at  the  65%  semispan  location  also  shows 
a  double  shock,  but  with  the  separation  distance  between  the  two  shocks 
decreasing.  Finally,  Figure  16.c  shows  the  pressure  distribution  at  95% 
semispan  location,  where  only  a  single  shock  exists  in  the  flow.  The 
configuration  obviously  results  in  the  lambda  double  shock  pattern  for 
transonic  conditions  on  a  swept  wing,  where  the  two  shocks  coalesce  to 
form  a  single  shock  near  the  wing  tip.  The  results  follow  the  same 
pattern  as  the  RAE  airfoil,  where  the  Roe  scheme  consistently  models 
shocks  in  much  fewer  points  with  less  smearing.  None  of  the  schemes  show 
any  sign  of  flow  separation,  which  also  agrees  with  the  experimental 
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data.  Figure  17  shows  a  series  of  pressure  contours  along  the  wing  span. 
The  contours  in  Figure  17. a  at  44%  semispan  location  clearly  show  a  weak 
shock  at  the  20%  chord  position  and  a  much  stronger  shock  near  the  60% 
chord  position.  Figure  17.b  shows  the  result  at  65%  semispan  location. 
Strong  shocks  exist  at  about  the  20%  and  50%  chord  positions,  and  the 
rear  most  shock  is  located  further  forward.  The  95%  semispan  location 
(Figure  17.c)  clearly  shows  the  shocks  having  coalesced  to  form  one  at 
the  25%  chord  position  and  this  shock  is  by  far  the  strongest  shock  of 
all  those  observed  in  Figure  17.  Figure  17.d  shows  a  view  of  the 
contours  along  the  upper  surface  and  the  double  shock  pattern  coalescing 
into  a  single  shock  at  the  tip.  These  solutions  compare  well  with  the 
experimental  data.  The  comparison  of  the  Roe  scheme  with  the  EMULE 
scheme  is  as  expected.  The  Roe  scheme  shows  slightly  improved  results  in 
the  leading  edge  expansion  region  and  directly  downstream.  The  Roe 
scheme  is  obviously  superior  again  in  all  shock  regions  where  BMULE  has 
smeared  the  shocks  over  several  grid  cells.  There  is  little  or  no 
difference  between  the  two  solutions  for  the  lower  surface  calculations. 
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SECTION  VI 


UNSTEADY  RESULTS 

6.1  Oscillating  NACA  0012  Airfoil 

Unsteady  two-dimensional  calculations  were  performed  for  the  NACA 
0012  airfoil  oscillating  in  pitch  about  the  25%  chord  point  as  described 
in  paragraph  4.4.  The  grid  used  was  the  same  as  previously  presented  in 
Figure  1.  The  calculations  presented  in  paragraph  4.4  were  all  for  only 
the  first  60  degrees  of  oscillatory  motion.  The  calculations  were  all 
stepped  at  this  point  in  the  cycle  to  compare  results  from  different  time 
step  sizes  and  subiteration  combinations.  The  results  to  be  presented 
here  are  exanples  of  taking  selected  cases  from  those  presented 
previously  and  continuing  the  motion  for  a  full  4  cycles  (1,440  degrees 
of  oscillatory  motion)  and  catparing  the  lift  and  moment  coefficient  time 
histories  and  selected  "snap-shots"  of  the  pressure  distributions  from 
the  last  cycle  of  motion. 

The  first  case  presented  is  a  comparison  of  4  minimum  time  step 
subiterations  and  4  local  time  step  subiterations  both  for  DT  =  0.05111, 
which  corresponds  to  a  maximum  CFL  of  45,000  compared  with  a  0 
subiteration  case  using  a  DT  =  0.005111  time  step.  The  0  subiteration 
case  at  DT  =  0.005111  was  used  at  the  small  time  step  for  comparison 
purposes  instead  of  the  DT  =  0.00102,  which  was  used  for  the  first  60:' 
motion  cases.  This  larger  time  step  was  selected  for  cost  savings,  since 
to  calculate  4  cycles  of  motion  with  the  small  time  step  (DT  =  0.00102) 
would  require  200,000  iterations,  while  the  larger  (DT  =  0.005111) 
required  40,000  iterations.  The  lift  coefficient  versus  time  plot  is 
shown  in  Figure  18. a  and  the  moment  versus  time  shown  in  Figure  18.b. 
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Only  a  slight  difference  in  the  lift  curve  is  noticeable  at  the  maximum 
and  minimum  values  of  lift  due  primarily  to  the  slight  displacement  of 
the  shock  location  resulting  form  the  larger  time  steps.  The  difference 
is  more  pronounced  for  the  moment  curves  (Figure  18. b) . 

The  second  case  presented  is  a  comparison  of  a  one-block  grid  and  a 
three-block  grid  shown  in  Figure  19  for  the  NACA  0012  airfoil.  Four 
minimum  time  step  subiterations  were  used  for  both  the  one-block  and  the 
three-block  grids  with  a  time  step  of  DT  =  0.05111.  Figure  20. a  -  20. h 
show  snap-shots  in  time  of  unsteady  pressure  distributions  for  the  two 
blocking  arrangements .  The  slight  differences  in  sere  of  the  pressure 
distributions  can  be  attributed  to  a  couple  of  reasons.  First, 
as  mentioned  in  Section  IV  the  boundary  conditions  used  at  block 
interfaces  were  those  developed  by  Be lk,  and  the  error  generated  by 

approximating  the  &Q  (0(  At'/^x) )  causes  a  slight  misplacement  of 
the  shock  position  between  the  one  and  three  block  cases,  especially  for 
the  times  when  the  shock  is  traveling  the  fastest  (i.e.,  Figures  20. a  and 
20.e).  Secondly,  the  use  of  blocks  degrades  the  convergence  rates  for 

steady  solutions  ( see  Belk  2 )  and  as  a  result  the  subiterations  for  the 
three-block  case  do  not  converge  as  rapidly  as  do  the  subiterations  for 
the  one-block  case,  which  accounts  for  the  slight  differences  in  the 
unsteady  pressure  distributions  elsewhere  in  the  solution. 

Figure  21  shows  the  pressure  contours  for  the  three-block  case  at 
the  end  of  the  fourth  cycle  of  motion.  Note  that  even  though  the  airfoil 
is  back  to  the  0  degree  angle  of  attack  position,  the  flow  field  is  not 
syrrmetric  and  a  relatively  strong  shock  still  exists  on  the  lower 
surface.  This  is  a  good  example  of  how  the  aerodynamics  lag  the 
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time-dependent  motion  of  the  body. 

The  third  and  last  two-dimensional  case  presented  is  a  comparison  of 
the  same  4  minimum  time  step  subiteration  case  using  the  first  order 
backward  Euler  differencing  formula,  Equation  (3.5a),  cat  pared  with  4 
minimum  time  step  2nd  order  subiteration  using  the  three -point  backward 
difference  formula  given  in  Equation  (3.5b).  Figure  22  show  snap-shots 
of  the  unsteady  pressures  at  the  same  times  as  shown  in  Figure  20, 
compared  with  the  experimental  data27 .  As  noted  by  Belk  2  this 
data  is  somewhat  suspect  due  to  obvious  asymmetric  properties,  even 
though  the  solution  should  be  very  nearly  a  symmetric  solution.  Belk 
determined  a  steady  angle  of  attack  somewhat  larger  than  the  0.016° 
should  be  used  to  improve  the  correlation  between  experimental  data  and 
calculations . 

The  second  order  solutions  show  the  shock  leading  the  shock 
location  predicted  by  the  first  order  solutions  for  ET  =  0.005111  and  DT 
=  0.05111  with  4  minimum  subiterations.  The  largest  difference  occurs 
when  the  shock  is  traveling  the  fastest,  such  as  the  70°  of  oscillatory 
motion  (Figure  22.b).  Note  that  at  the  point  where  the  shock  is 
basically  stationary  (i.e.,  160°  of  motion)  the  two  CT  =  0.05111  cases 
for  first  and  second  order  schemes  are  essentially  identical.  At  all  the 
times  shown  in  Figure  22  the  second  order  solution  shows  a  closer 
agreement  with  the  experimental  data  than  does  the  two  first  order 
solutions . 

The  comparison  of  the  2nd  order  solution  and  the  experimental  data 
indicates  an  even  closer  agreement  between  the  two  would  be  possible  if  a 
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slightly  larger  steady  angle  of  attack  was  used  for  the  calculations. 

For  exarrple,  a  larger  steady  angle  of  attack  added  to  the  unsteady  angle 
of  attack  would  in  Figures  22. b  -  22 .d  result  in  a  larger  expansion 
region  and  a  shock  location  farther  aft  as  indicated  by  the  experimental 
data.  Likewise,  in  Figures  22. f  -  22. h  when  the  airfoil  is  at  negative 
angle  of  attack  a  larger  (more  positive)  steady  angle  of  attack  would 
result  in  a  smaller  expansion  region  on  the  lower  surface  and  the  shock 
located  farther  forward  as  indicated  by  the  experimental  data.  This 
confirms  Be Ik's  findings  on  the  same  subject. 

6.2  Oscillating  Supercritical  Rectangular  Planform  Wing 

The  supercritical  rectangular  planform  wing  was  used  as  the  three- 
dimensional  unsteady  check  case  for  this  work  due  to  the  extensive  amount 
of  wind  tunnel  test  data35  that  exist.  The  wing  has  an  aspect  ratio 
of  2  and  was  oscillated  about  the  wing  pitch  axis  located  at  the  40% 
chord.  The  experimental  data  used  was  taken  in  Freon,  and  therefore  the 
ratio  of  specific  heats  is  y  =  1.131  with  a  transonic  Mach  number  of 
0.7  and  a  steady  angle  of  attack  of  4°. 


The  grid  for  this  case  had  221  x  40  x  15  points.  The  grid  used  a 
"C"  mesh  in  the  streamwise  direction  and  a  "H"  mesh  in  the  spanwise 
direction.  Each  of  the  K  planes  (constant  spanwise  location)  have  the 
same  point  distribution  as  the  two-dimensional  grid  used  for  the  NACA 
0012  airfoil  calculations,  Figure  1 .a,  with  clustering  at  the  leading  and 
trailing  edges.  The  distribution  was  selected  to  give  the  most 
orthogonal  grid  possible.  Figure  2 3. a.  There  are  10  such  planes  of  data 
distributed  down  the  wing  span,  with  the  planes  being  clustered  near  the 
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tip.  The  wing  tip  was  not  modeled  in  this  case  as  was  done  for  the  ONERA 
M6  wing,  but  instead  the  airfoil  was  sirrply  collapsed  to  a  line  at  the 
K= 1 1  plane,  which  was  placed  just  out  from  the  wing  tip.  An  additional  4 
grid  planes  were  distributed  in  the  spanwise  direction  to  the  far field 
boundary  (Figure  23. b).  The  spacing  for  the  first  point  off  the  wing  was 
such  that  a  minimum  ¥  of  1.5  occurred  on  the  wing. 

The  grid  described  was  used  for  steady  calculations  to  compare  with 
the  experimental  data  from  Reference  (35)  to  verify  the  quality  of  the 
grid  before  performing  the  unsteady  calculations.  The  results  near  the 
wing  tip  cure  of  particular  interest  since  the  grid  used  did  not  model  the 
tip  correctly.  The  transonic  flow  conditions  were  Mach  =  0.701  and  4.0° 
angle  of  attack,  which  was  considerably  different  from  the  supercritical 
design  conditions  of  Mach  =  0.8  and  0.0°  angle  of  attack,  and  results 
in  a  sharp  shock  near  the  20%  chord.  The  Figures  24. a  -  24. d  show 
excellent  agreement  with  the  experimental  data  and  confirm  the  quality  of 
the  grid  (including  the  collapsed  wing  tip)  to  accurately  model  the  flow 

field  for  the  rectangular  wing.  Be  Ik  2  has  reported  a  predicted  shock 
position  downstream  of  the  position  shown  in  the  experimental  data  using 
an  inviscid  approximation.  The  inclusion  of  the  viscous  terms  has 
corrected  the  position  of  the  predicted  shock  to  more  nearly  coincide 
with  the  experimental  data. 

The  three-dimensional  unsteady  calculations  were  obtained  at  a  Mach 
number  0.699  and  a  steady  angle  of  attack  of  4.03°.  The  angle  of 
attack  was  varied  according  to  Equation  (4.20)  with  an  unsteady  anplitude 
of  1.03S5  and  a  reduced  frequency.  Equation  (4.19),  of  0.358. 
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A  steady  starting  condition  was  obtained  using  a  CFL  of  5  for  1 , 000 
iterations,  after  which  the  notion  was  impulsively  started.  The 
calculations  were  run  for  three  cycles  of  notion.  During  the  third  cycle 
the  pressure  coefficients  were  saved  at  each  time  step.  The  time  step 
size  was  DT  =  0.069,  which  corresponds  to  360  time  steps  per  cycle  of 

notion  (same  as  used  by  Be lk 2  for  the  inviscid  calculations),  and  5 
minimum  time  step  subiterations  were  used  to  converge  the  solution  at 
each  time  step  by  nearly  one  order  of  magnitude.  The  maximum  CFL  for 
this  case  was  approximately  17,000. 

The  magnitude  and  phase  of  the  unsteady  pressure  coefficients  from 
the  third  cycle  of  notion  was  obtained  by  Fourier  analysis.  The 
magnitude  and  phase  of  the  unsteady  pressure  coefficients  are  presented 

along  with  the  inviscid  results  from  Belk  2  and  the  experimental 
data35,  for  the  semispan  locations  in  Figures  2 5. a  -  25.g.  Again 

the  predicted  results  compare  well  with  the  experimental  data.  Belk  2 
reported  a  misplacement  of  the  sharp  spike  in  magnitude  of  the  unsteady 
pressure  coefficient  due  to  the  inviscid  code  misplacing  the  shock.  Once 
again  the  inclusion  of  the  viscous  terms  has  seemingly  corrected  this 
problem  to  give  good  agreement  with  the  experimental  data. 

The  good  comparisons  tend  to  drop  off  as  the  wing  tip  is  approached, 
such  that  at  the  95%  semispan  location  the  phase  comparison  with  the 
experimental  data  (Figure  25. f  and  25. g)  is  not  as  good.  The  shift  in 
the  phase  for  the  upper  surface  is  obviously  misplaced,  and  the  phase  for 
the  lower  surface  is  consistently  under  predicted.  The  story  is  similar 
but  not  as  severe  at  85%  semispan  location.  It  should  be  noted  that  the 
experimental  data  at  this  station  for  the  lower  surface  is  highly 
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suspect,  since  the  data  does  not  follow  the  obvious  trend  of  the  data  at 
the  stations  on  either  side  (inboard  or  outboard).  The  reduced  accuracy 
at  the  tip  should  ccme  as  no  great  surprise  considering  the  grid 
treatment  of  the  wing  tip.  Be Ik's  2  results,  while  showing  some  decline 
in  accuracy  near  the  tip,  do  not  show  as  dramatic  a  change.  This 
strengthens  the  argument  for  blaming  the  grid  topology  used  ("C-H")  since 
Be lk  used  a  ("C-0")  type  grid  similar  to  the  one  used  in  this  work  for 
the  ONERA  wing. 
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SECTION  VII 


CONCLUSIONS 

Two  algorithms  have  been  presented  for  corputing  unsteady  thin-layer 
Navier -Stokes  solutions  for  oscillating  airfoils  and  wings.  The 
algorithms,  flux  vector  split  (FVS)  and  flux  difference  split  (FD6)  were 
first  compared  for  steady  viscous  calculations.  The  computed  laminar 
boundary  layer  profiles  for  a  flat  plate  were  compared  with  the  Blasius 
solutions  and  the  FDS  was  the  clear  winner  since  it  managed  to  capture 
the  laminar  profile  with  only  three  grid  points  inside  the  boundary 
layer,  while  with  four  times  as  many  grid  points  inside  the  boundary 
layer,  the  FVS  scheme  could  still  not  match  the  comparison  of  FDS  with 
the  Blasius  solutions.  The  inability  of  the  FVS  scheme  to  capture  a 
laminar  boundary  layer  is  due  to  the  excessive  numerical  viscosity  in  the 
scheme.  Similar  improvements  in  the  computed  steady  solutions  for  the 
RAE  2822  airfoil  and  the  ONERA  M6  wing  were  observed  for  the  FDS  scheme. 
The  numerical  dissipation  in  the  FVS  scheme  is  also  observed  in  shocks  in 
the  form  of  smearing  the  shock  waves  over  several  grid  cells.  The 
improved  solutions  for  the  FDS  scheme  did  not  come  cheap,  since  the  FDS 
scheme  has  a  higher  operational  count  than  does  the  FVS  (approximately 
20%  higher),  and  the  convergence  rate  for  the  FDS  scheme  is  slower  than 
the  FVS  scheme,  therefore  requiring  more  iterations  to  reach  the  same 
convergence  level  (approximately  25%  more  iterations). 

The  second  order  FES  scheme  with  the  minmod  limiter  was  selected  as 
the  algorithm  to  be  used  for  the  unsteady  thin-layer  Navier -Stokes 
computations.  The  emphasis  was  placed  on  obtaining  time-accurate  viscous 
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solutions  for  oscillating  airfoils  and  wings  as  cheaply  as  possible. 

Since  the  diffusive  terms  were  treated  explicitly,  the  formal  temporal 
accuracy  of  the  algorithm  without  subiterations  was  limited  to  1st  order 
accurate.  A  form  of  Newton  subiterations  was  implemented  to  converge  the 
unsteady  calculations  at  each  time  step  before  progressing  the  solution 
to  the  next  time  step.  The  use  of  these  subiterations  then  allowed  the 
calculations  to  progress  at  much  larger  time  step  sizes  (and  CFL  number) 
than  would  otherwise  be  possible.  The  use  of  subiterations  also  provided 
the  capability  for  a  2nd  order  time  accurate  scheme  at  convergence  of  the 
subiterations . 

The  use  of  subiterations  were  shown  to  significantly  improve  the 
quality  of  the  solution,  and  in  seme  cases  allow  a  solution  to  be 
obtained  which  otherwise  would  have  been  impossible  at  the  time  step  size 
being  used.  Unsteady  calculations  for  the  two-dimensional  NACA  0012 
airfoil  were  performed  to  help  select  the  proper  combinations  of  the  step 
size  and  subiteration  number  to  give  the  best  possible  solution.  The 
results  from  these  calculations  were  used  to  make  more  extensive  two  and 
three-dimensional  unsteady  calculations  to  compare  with  experimental 
data.  The  three-dimensional  calculations  for  the  supercritical 
rectangular  planform  wing  showed  excellent  agreement  with  the 
experimental  data,  and  the  use  of  subiterations  with  the  FIDS  scheme 
resulted  in  a  relatively  efficient  algorithm  which  could  be  used  in  more 
complex  three-dimensional  unsteady  viscous  calculations. 
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TABLE  I 


SECOND  ORDER  TRUNCATION  ERROR 


4* 

NAME 

b 

2nd  order  TE 

1/3 

Third-Order 

4 

0 

-1 

Fully  Upwind 

2 

^(ax)2 

0 

Frcrm's 

3 

^(Ax)2 

1/2 

Low  d  order 

5 

-^(AX)2 

1 

Central 

00 

-g(Ax)2 

-1/3 

No  Name 

5 

1 

^(ax)2 

TABLE  II 


SUMMARY  OF  TIME  STEP  SIZE  RUNS  ' 


TIME  STEP 

#  ITERATIONS 

CPU  SECONDS 

#  SUBS 

0.0010 

8250 

46060 

0 

0.0051 

1640 

9161 

0 

0.0102 

833 

737 

0 

0.0511 

167 

149 

0 

0.1022 

83 

76 

0 

U.2044 

42 

39 

0 

*  First  60P  of  motion. 
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TABLE  III 


SUMMARY  OF  SUBITERATION  RUNS  ' 


TIME  STEP 

#  ITERATIONS 

CPU  SECONDS 

#  SUBS 

0.0511 

167 

1546 

16  MIN 

0.0511 

167 

497 

4  MIN 

0.0511 

167 

324 

2  MIN 

0.2044 

42 

729 

32  MIN 

0.2044 

42 

383 

16  MIN 

0.2044 

42 

214 

8  MIN 

0.0511 

167 

1530 

-- 

16  LOCAL 

0.0511 

167 

495 

4  LOCAL 

0.0511 

167 

322 

2  LOCAL 

0.2044 

42 

735 

32  LOCAL 

0.2044 

42 

382 

16  LOCAL 

0.2044 

42 

208 

8  LOCAL 

’  First  60°  of  motion. 
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FIGURES 


70 


Close-up  of  NACA  0012  section 


NACA  0012  UNSTEADY  PRESSURE  DISTRIBUTION 

M A  CH  =  Q  .  755  RE^.  5E06  ALPHA  =  0.016 
(ZERO  SUB ! TERAT I ONS ) 


Unsteady  pressure  distribution  for  different 


Unsteady  pressure  distributions  using  minimum  time  step 
subiterations  for  first  60  degrees  of  motion  DT  =  .0511 


Unsteady  pressure  distributions  using  minimum  time  step 
subiterations  for  first  60  degrees  of  motion  DT  =  .20^ 


NACA  0012  UNSTEADY  PRESSURE  DISTRIBUTION 

M  A  C 11  =  0  .  7  5  5  KE  =  5.5E0G  A  L  P 11 A  =  0  .  0  L  U 
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Figure  5.  Unsteady  pressure  distributions  using  local  time  step 
subiterations  for  first  60  degrees  of  motion  DT  =  .0511 


Unsteady  pressure  distributions  using  both  min  and  local 
subiterations  for  first  60  degrees  of  motion  DT  =  .0511 


Unsteady  pressure  distributions  using  1st  and  2nd  order 
schemes,  16  min  time  step  subiterations 


LAMINAR  FLAT  PLATE  VELOCITY  DISTRIBUTION 

M A C H  =  0  .  5  RE L=  i  0 , OOO 


Figure  9.  Velocity  profiles  for  laminar  flat  plate 


E  2822 


2822  STEADY  PRESSURE  D1 

C  A  S  E  =  9  M  A  C  U  =  0  .  7  \i  4  K  K  =  0  .  5  E  0  0  A 
(HOE  -  M1NMOD  2rd  OKDEU)  (UMULE  - 


8  2  2  BOUNDARY  LAYER  PROF!  LE 

9  M  A  C  H  =  0 . 7  3  4  RE=G.5E06  ALPHA=2.79 

MINMOD  2nd  ORDER)  (BMULE  -  2nd  ORDER) 


Boundary  layer  velocity  profiles  at  32$  chord  (upper  surface) 


822  SKIN  FRICTION 

3  =  0  .  7  34  KE-6  .  5E06  ALPHAS. 79 


surface  skin  friction 


R A E  2  8  2  2  BOUNDARY  LAYER  PROF 


Boundary  layer  velocity  profiles  at  32$  chord  (upper  surface) 


RAE  2822  BOUNDARY  LAYER  PROFILE 

C  A  S  E  =  9  MAC  II  =  0  .  734  RE  =  0.5EQ(3  A  L I'  li  A  =  2  .  7  9 
(ROE  -  MINMOD  3rd  ORDER)  (13MULE  -  2nd  ORDER 


Boundary  layer  velocity  profiles  at  57$  chord 


RAE  2022  SKIN  FRICTION 

CASE=9  MACH=0 . 734  RL'  =  0.5E0C  ALPHA=2.79 
ROE  -  M I NMOD  3rd  ORDER)  (DMULE  -  2ud  ORDER) 


surface  skin  friction 


2822  STEADY  PRESSURE  D  I  ST R 

CASE  =  9  MACH  =  0.734  RE  =  (i.5E0C  A  L  P  U  A  = 
ROE  -  SUPDEE  2rd  ORDER)  (DMULE  -  2nd 


RAE  2822  BOUNDARY'  LAYER  PROFILE 

C  A  S  E  =  9  MA  Cfi  =  0  .  734  RE=G.5EOC  ALPHA=2.79 


Figure  13.  (continued) 


822  SKIN  FRICTION 

Si=0.7  34  RE  =  6  .  5EQ6  ALPHA=2.79 
2nd  ORDER)  ( BMULE  -  2nd  ORD1 


surface  skin  friction 
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Figure  14.  Convergence  histories  for  Roe  and  BMULE  schemes 


Figure  15.  (concluded) 


ONERA  M- 6  WING  STEADY  PRESSURE  DISTRIBUTION 

MA CH=  O . 8 4  RE  =  2 . GEOG  ALPHA  =  3.0G  WING  STATION  =  6 

(  BMULE  vs  ROE  ) 


semispan  location 


Upper  surface  pressure  contours 


LIFT  COEFFICIENT  NACA  0012 

MACH=0 .755  MEAN  ALPHA=0.0I6  DALPHA=2.51  K=O.1620 


105 


Figure  18.  Comparisons  of  lift  and  moment  coefficients  using 
4  min  and  4  local  subiterations 


MOMENT  COEFFICIENT  NACA  0012 

M  ACH  =  0 . 755  MEAN  ALPHA=0.016  DALPHA=2.5L  K- 0  .  16  2  8 
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Moment  coefficient  vs  time 
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Figure  19.  Three  block  version  of  "C"  grid  for  NACA  0012 


NACA  0012  UNSTEADY  PRESSURE  DISTRIBUTION 

MACn=0 . 755  R  E  =  5 . 5E06  ALPHA=0.0I6  K=0.JC20 


Unsteady  pressure  distributions  for  1  and  3  blocks 


12  UNSTEADY  PRESSURE  DISTRIBUTION 

C  [I  =  O  .  7  5  5  R  E  =  5  .  5E06  ALPHA=0.016  K  =  0.1620 


degrees  of  oscillatory  motion,  AOA  increasing 


PRESSURE 


oscillatory  motion,  AOA  decreasing 


N A C A  0012  UNSTEADY  PRESSURE  DISTRIBUTION 

M  A  C  U  =  0 . 755  R  E  =  S . 5  E  0  6  ALPUA=0.01G  K  =  0  .  i G  2  0 


Figure  20.  (continued) 


NACA  0012  UNSTEADY  PRESSURE  DISTRIBUTION 

UACn=0 . 755  R  E  =  5 . 5  E  0  6  ALPHA=0.016  K=0.1G28 


degrees  of  oscillatory  motion,  AOA  decreasing 


UNSTEADY  PRE.SSURE  DISTRIBUTION 


degrees  of  oscillatory  motion,  AOA  increasing 


NACA  0012  UNSTEADY  PRE-SSURE  DISTRIBUTION 

MACH  =  0.755  RE=5.5E06  ALPBA=0.016  K  =  0 . 1 6  2  8 
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Figure  22.  Unsteady  pressure  distributions  for  1st  and  2nd  order 
time  accuracy  using  H  min  subiterations  and  DT  =  .0511 


STEADY  PRESSURE  DISTRIBUTION 

3  RE  =  5 . 5E0C  ALP1IA  =  0.016  K  =  0.1CE0 


12  UNSTEADY  PRESSURE  DISTRIBUTION 

CH  =  O  .  7  6  5  R  E  =  5  .  5  £  D  6  ALPIIA  =  0.016  K  =  0.1628 


degrees  of  oscillatory  motion 


NACA  0012  UNSTEADY  PRESSURE  DISTRIBUTION 

MACH=Q .755  RE=5.5E0G  ALPUA.  =  0  .  0  1  6  K=0.1628 


degrees  of  oscillatory  motion 


UNSTEADY  PRESSURE  DISTRIB 
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degrees  of  oscillatory  motion 


NACA  0012  UNSTEADY  PRESSURE  DISTRIBUTION 


Figure  22.  (concluded) 


RECTANGULAR  SUPERCRITICAL  WING 

STEADY  PRESSURE  DISTRIBUTION 


Steady  pressure  distributions  for  the  rectangular 
supercritical  wing  compared  with  experiment 


CTANGULAR  SUPERCR 

STEADY  PRESSURE  D  I  S‘I 


semispan  calculated,  58%  experiment 


CTANGULAR  SUPERC 

STEADY  PRESSURE  DI 


semispan  calculated,  8 2%  experiment 


RECTANGULAR  SUPERCRITICAL  WING 

STEADY  PRESSURE  DISTRIBUTION 
MACQ=0 . 7  0  4  REM.03EQ6  ALPBA  =  4.00 


aemiapan  calculated,  95H  experiment 


RECTANGULAR  SUPERCRITICAL  WING 

UNSTEADY  PRESSURE  DISTRIBUTIONS 
MACH=0 . 699 ;  4  DEG  STEADY  ALPHA;  1  DEG  UNSTEADY  ALPHA,  k=0.35G 
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Magnitude  at  3^$  semispan 


RECTANGULAR  SUPERCRITICAL  WING 

UNSTEADY  PRESSURE  DISTRIBUTIONS 
.699;  4  DEG  STEADY  ALPHA;  1  DEG  UNSTEADY  ALPHA,  k=0.S58 


continued) 


CTANGULAR  SUPERCRITICAL  WING 

UNSTEADY  PRESSURE  DISTRIBUTIONS 
l;  4  DEG  STEADY  ALPHA;  1  DEG  UNSTEADY  ALPUA,  k=0 . 35B 


Magnitude  at  61$  semispan 


RECTANGULAR  SUPERCRITICAL  WING 

UNSTEADY  PRESSURE  DISTRIBUTIONS 
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(continued ) 


RECTANGULAR  SUPERCRITICAL  WING 

UNSTEADY  PRESSURE  DISTRIBUTIONS 
.699;  4  DEG  STEADY  ALPHA;  1  DEG  UNSTEADY  ALPHA,  k=0 .358 
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RECTANGULAR  SUPERCRITICAL  WING 

UNSTEADY  PRESSURE  DISTRIBUTIONS 
MACQ=0 . 699  ;  4  DEG  STEADY  ALPHA;  1  DEG  UNSTEADY  ALPHA,  k  =  0.35U 


(continued) 


CTANGULAR  SUPERCRITICAL  WING 

UNSTEADY  PRESSURE  DISTRIBUTIONS 
i;  4  DEG  STEADY  ALPHA;  1  DEG  UNSTEADY  ALPHA, 


(continued) 


RECTANGULAR  SUPERCRITICAL  WING 

UNSTEADY  PRESSURE  DISTRIBUTIONS 
MACU=0 . 699 ;  4  DEG  STEADY  ALPHA;  1  DEG  UNSTEADY  ALPHA,  lc=0.358 


Figure  25.  (concluded) 


APPENDIX  A 


CURVILINEAR  TRANSFORMATION 

The  Navier-Stokes  equations  in  Cartesian  coordinates  from  Elquation 
(2.7)  are 


3?  +  3f  +  ag  +  =  o 

0t  0x  0y  0z 

The  general  curvilinear  coordinates  are  given  as 

£  =  (x,y,z,t ) 

1]  =  i]  (x,y,z,t ) 

£  =  £  (x,y,z,t ) 

T  =  T(t) 


(A.  1  ) 


(A. 2) 


Expanding  the  Cartesian  coordinates  in  terms  of  the  curvilinear 
coordinates  using  the  chain  rule  gives 


0  8t  3  +  03;  0  +  01]  8  +  3C  9 

0t  0t  0T  0t  0H  at  0)-j  0t  0C 

<L  =  ;K  <L  +  ;h  ^  +K  <L 

0X  0x  0t-  0x  0,1  0x  0^ 

±_  =  §£  a_  +  an  a_  +  3C  a_ 

0y  0y  0£  0y  0i)  0y  a<r 

a_  =  aia_  +  a]ijL+KiL 

02  02  0£  0z  0,i  0 z  0C 


The  inverse  of  this  relationship  is  determined  by  expanding  the 
curvilinear  coordinates  in  terms  of  Cartesian  coordinates  (in  matrix 
form)  as 
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V  z 

yT  T 


y?  z5 


V  z„ 

•  11  11 


z 


c  c 


(A. 4) 


Let  Equation  (A.  4)  be  written  as 


[  ]  =  [A]"1  [  ] 


The  jacobian  of  the  transformation  is  given  by  the  det  of  [A]  defined  as 


J’  =  |A| 


z?y, ,XC  ■  ySxnzC  ‘  Wf] 


tj  J 


where, 


J  =  tx?Vr '  Vc'*VVc  "Vr)  +  z?Vr  "Vr,] 


( A- 5 ) 


The  metric  coefficients  are  determined  by  inverting  the  martix  A  to  yield 


v  r  VrVf’ 

^  -  j"  (ZnXr  ■ 

h  'J"  Vc  'y„V 


=  J'1  (Z^  -  > 

=  J  1  (X,-Zr  -  Z>-Xr  > 

s  V.  s  V 

=  r'  (ycxr  -  xpyr» 
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yj"  y>i~¥V 

y J"  VirVn1 

C,  =  J’1  (xvy  -  yex  > 

Z  Cl]  y 

T‘y 

*t  -Tt  ('*T«i_Vy  'Mz* 

"t  =Tt  XT1|x  "  ~  Zt'*Z  * 

r,  =  t,  <-XTCx-yTry-zT?z> 


Substituting  the  expanded  derivatives  in  Cartesian  coordinates  (Equation 
(A- 3 ) )  into  the  Navier-Stokes  equations  (Equation  (A. 1 ) )  results  in 


j  <!2  +  c  3a  +  „  ^a  +  (-  5a  + 

Tt  3t  st  a?  'H  an  m  ac 

c  af  . ,,  ar  .  r  ar  . 

?x  a?  "x  ^  +  rx  %  * 

sy  a?  V  an  ^y  ar 

?Z  af  +  "z  ali +  rz  af  =  0 


(A. 7) 


Then  by  multiplying  Equation  (A. 7)  by  J'  and  using  the  chain  rule 
Equation  (A.1)  becomes 


3Q  +  3F  .  3G  .  3H 
3t  3^  3 ii  3C  " 


(A. 8) 


where. 
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U  = £  +Hu+£v+£w 

^  t  ^x  ^y  ^z 

V  =  nt  +  \n  +  Hyv  +  H  w 

w  =  r  +(u  +  cv  +  Cw 

5  t  s  x  s  y  5  z 

Tf,  =  5J„  +  5yTyl 

T  =?T  +F  T  +?T 

?y  sx  xy  sy  yy  sz  xy 


T„,  =  Vx  +  'V1,.  +  V* 
T„  *  V.y  +  V»v  +  V-y 

Tnz  '  V.z  +  ’l,T,z  +  Vzz 

Tcz  =?,Tz<  +f,Ty»  +  tT.z 
Tt,  ‘  f.Tz,  +  f»Tyy  +  fzTz, 

Tfz  "  Cj„  +  fyTyz  +  ?zTzz 

rF  '  ?  q  +  E  q  +  ?  q 

s  ^x  *x  -y^y  ^z  z 

r„  "  \%  +  I)yqy  +  n,qz 
rr  -  r,q,  +  cyq,  +  C:q: 


(A. 10) 


The  viscous  stress  tensors  for  the  full  Navier -Stokes  equations  are 
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7xx  =  -  Vs+V^yV  “  ^‘VW'c  )] 

7yy  =  TT^  [2(?yV?+T,yVn4i:yVr  }  “  {*A'h'AAA)  '  (*z  VV'n'V'r  ^ 

7zz  =  [2(?2wp+\wn^zwr }  -  ^VVV^V  -  ^VV^yV3 

Txy  -  ^  ^  ♦y’n  +  fy^r  +  ^xve  +  \vn  +  rxvr] 

Txz  =  X  ^zS  +  ^zUr  +  ?::w5  +  T'XW-|  + 

jj  M 

Tyz  =  if  C^V?  +Vn  +  CrVr  +  ?yW?  +  \Wn  +  VV ] 

and  the  heat  flux  terms  are 

<1,  -  ^  Me  +  ,ST,  +  r,Tr] 

%  ■  y^T  Me  +  "yT'l  +  CyTf] 

<4  -  ^  Me  +  M,  +  t  Tr] 
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